1function [o,count,SSQ] = sumskipnan(x, DIM, W)
2% SUMSKIPNAN adds all non-NaN values.
3%
4% All NaN's are skipped; NaN's are considered as missing values.
5% SUMSKIPNAN of NaN's only  gives O; and the number of valid elements is return.
6% SUMSKIPNAN is also the elementary function for calculating
7% various statistics (e.g. MEAN, STD, VAR, RMS, MEANSQ, SKEWNESS,
8% KURTOSIS, MOMENT, STATISTIC etc.) from data with missing values.
9% SUMSKIPNAN implements the DIMENSION-argument for data with missing values.
10% Also the second output argument return the number of valid elements (not NaNs)
11%
12% Y = sumskipnan(x [,DIM])
13% [Y,N,SSQ] = sumskipnan(x [,DIM])
14% [...] = sumskipnan(x, DIM, W)
15%
16% x	input data
17% DIM	dimension (default: [])
18%	empty DIM sets DIM to first non singleton dimension
19% W	weight vector for weighted sum, numel(W) must fit size(x,DIM)
20% Y	resulting sum
21% N	number of valid (not missing) elements
22% SSQ	sum of squares
23%
24% the function FLAG_NANS_OCCURED() returns whether any value in x
25%  is a not-a-number (NaN)
26%
27% features:
28% - can deal with NaN's (missing values)
29% - implements dimension argument.
30% - computes weighted sum
31% - compatible with Matlab and Octave
32%
33% see also: FLAG_NANS_OCCURED, SUM, NANSUM, MEAN, STD, VAR, RMS, MEANSQ,
34%      SSQ, MOMENT, SKEWNESS, KURTOSIS, SEM
35
36
37%    This program is free software; you can redistribute it and/or modify
38%    it under the terms of the GNU General Public License as published by
39%    the Free Software Foundation; either version 3 of the License, or
40%    (at your option) any later version.
41%
42%    This program is distributed in the hope that it will be useful,
43%    but WITHOUT ANY WARRANTY; without even the implied warranty of
44%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
45%    GNU General Public License for more details.
46%
47%    You should have received a copy of the GNU General Public License
48%    along with this program; If not, see <http://www.gnu.org/licenses/>.
49
50%   Copyright (C) 2000-2005,2009,2011 by Alois Schloegl <alois.schloegl@gmail.com>
51%   This function is part of the NaN-toolbox
52%   http://pub.ist.ac.at/~schloegl/matlab/NaN/
53
54
55global FLAG_NANS_OCCURED;
56
57if nargin<2,
58        DIM = [];
59end;
60if nargin<3,
61        W = [];
62end;
63
64% an efficient implementation in C of the following lines
65% could significantly increase performance
66% only one loop and only one check for isnan is needed
67% An MEX-Implementation is available in sumskipnan.cpp
68%
69% Outline of the algorithm:
70% for { k=1,o=0,count=0; k++; k<N}
71% 	if ~isnan(i(k))
72% 	{ 	o     += x(k);
73%               count += 1;
74%		tmp    = x(k)*x(k)
75%		o2    += tmp;
76%		o3    += tmp.*tmp;
77%       };
78
79if isempty(DIM),
80        DIM = find(size(x)>1,1);
81        if isempty(DIM), DIM = 1; end;
82end
83if (DIM<1), DIM = 1; end; %% Hack, because min([])=0 for FreeMat v3.5
84
85%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
86% non-float data
87%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
88if  (isempty(W) && (~(isa(x,'float') || isa(x,'double')))) % || ~flag_implicit_skip_nan(), %%% skip always NaN's
89	if ~isempty(W)
90		error('SUMSKIPNAN: weighted sum of integers not supported, yet');
91	end;
92	x = double(x);
93	o = sum(x,DIM);
94	if nargout>1
95		sz = size(x);
96		N  = sz(DIM);
97		sz(DIM) = 1;
98		count = repmat(N,sz);
99		if nargout>2
100			x = x.*x;
101			SSQ = sum(x,DIM);
102		end;
103	end;
104	return;
105end;
106
107if ~isempty(W) && (size(x,DIM)~=numel(W))
108	error('SUMSKIPNAN: size of weight vector does not match size(x,DIM)');
109end;
110
111%% mex and oct files expect double
112x = double(x);
113
114%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
115% use Matlab-MEX function when available
116%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
117%if 1,
118try
119
120	%% using sumskipnan_mex.mex
121	if issparse(x),
122		fprintf(2,'sumskipnan: sparse matrix converted to full matrix\n');
123		x = full(x);
124	end;
125
126	%% !!! hack: FLAG_NANS_OCCURED is an output argument, reserve memory !!!
127	if isempty(FLAG_NANS_OCCURED),
128		FLAG_NANS_OCCURED = logical(0);  % default value
129	end;
130
131	if (nargout<2),
132		o = sumskipnan_mex(real(x),DIM,FLAG_NANS_OCCURED,W);
133		if (~isreal(x))
134			io = sumskipnan_mex(imag(x),DIM,FLAG_NANS_OCCURED,W);
135			o  = o + i*io;
136		end;
137		return;
138	elseif (nargout==2),
139		[o,count] = sumskipnan_mex(real(x),DIM,FLAG_NANS_OCCURED,W);
140		if (~isreal(x))
141			[io,icount] = sumskipnan_mex(imag(x),DIM,FLAG_NANS_OCCURED,W);
142			if any(count(:)-icount(:))
143				error('Number of NaNs differ for REAL and IMAG part');
144			else
145				o  = o+i*io;
146			end;
147		end;
148		return;
149	elseif (nargout>=3),
150		[o,count,SSQ] = sumskipnan_mex(real(x),DIM,FLAG_NANS_OCCURED,W);
151		if (~isreal(x))
152			[io,icount,iSSQ] = sumskipnan_mex(imag(x),DIM,FLAG_NANS_OCCURED,W);
153			if any(count(:)-icount(:))
154				error('Number of NaNs differ for REAL and IMAG part');
155			else
156				o  = o+i*io;
157				SSQ = SSQ+iSSQ;
158			end;
159		end;
160		return;
161	end;
162end;
163
164if ~isempty(W)
165	error('weighted sumskipnan requires sumskipnan_mex');
166end;
167
168%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
169% count non-NaN's
170%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
171if nargout>1,
172        count = sum(x==x,DIM);
173	FLAG_NANS_OCCURED = any(count(:)<size(x,DIM));
174end;
175
176%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
177% replace NaN's with zero
178%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
179x(x~=x) = 0;
180o = sum(x,DIM);
181
182if nargout>2,
183        x = real(x).^2 + imag(x).^2;
184        SSQ = sum(x,DIM);
185end;
186
187%!assert(sumskipnan([1,2],1),[1,2])
188%!assert(sumskipnan([1,NaN],2),1)
189%!assert(sumskipnan([1,NaN],2),1)
190%!assert(sumskipnan([nan,1,4,5]),10)
191%!assert(sumskipnan([nan,1,4,5]',1,[3;2;1;0]),6)
192