Mercurial > forge

annotate extra/NaN/src/sumskipnan_mex.cpp @ 7889:c101c486d80a octave-forge

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
fix web address
author schloegl
date Thu, 27 Jan 2011 17:10:36 +0000
parents b9f35668b55e
children db5092052107
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
1 //-------------------------------------------------------------------
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
2 // C-MEX implementation of SUMSKIPNAN - this function is part of the NaN-toolbox.
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
3 //
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
4 //
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
5 // This program is free software; you can redistribute it and/or modify
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
6 // it under the terms of the GNU General Public License as published by
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
7 // the Free Software Foundation; either version 3 of the License, or
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
8 // (at your option) any later version.
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
9 //
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
10 // This program is distributed in the hope that it will be useful,
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
11 // but WITHOUT ANY WARRANTY; without even the implied warranty of
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
12 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
13 // GNU General Public License for more details.
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
14 //
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
15 // You should have received a copy of the GNU General Public License
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
16 // along with this program; if not, see <http://www.gnu.org/licenses/>.
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
17 //
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
18 //
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
19 // sumskipnan: sums all non-NaN values
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
20 // usage:
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
21 // [o,count,SSQ] = sumskipnan_mex(x,DIM,flag,W);
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
22 //
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
23 // SUMSKIPNAN uses two techniques to reduce errors:
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
24 // 1) long double (80bit) instead of 64-bit double is used internally
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
25 // 2) The Kahan Summation formula is used to reduce the error margin from N*eps to 2*eps
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
26 // The latter is only implemented in case of stride=1 (column vectors only, summation along 1st dimension).
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
27 //
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
28 // Input:
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
29 // - x data array
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
30 // - DIM (optional) dimension to sum
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
31 // - flag (optional) is actually an output argument telling whether some NaN was observed
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
32 // - W (optional) weight vector to compute weighted sum (default 1)
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
33 //
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
34 // Output:
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
35 // - o (weighted) sum along dimension DIM
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
36 // - count of valid elements
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
37 // - sums of squares
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
38 //
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
39 //
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
40 // $Id$
7889
c101c486d80a fix web address
schloegl
parents: 7888
diff changeset
41 // Copyright (C) 2009,2010,2011 Alois Schloegl <a.schloegl@ieee.org>
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
42 // This function is part of the NaN-toolbox
7889
c101c486d80a fix web address
schloegl
parents: 7888
diff changeset
43 // http://pub.ist.ac.at/~schloegl/matlab/NaN/
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
44 //
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
45 //-------------------------------------------------------------------
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
46
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
47
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
48
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
49
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
50 #include <math.h>
7888
b9f35668b55e replace <inttypes.h> with <stdint.h>
schloegl
parents: 7301
diff changeset
51 #include <stdint.h>
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
52 #include "mex.h"
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
53
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
54 inline int __sumskipnan2w__(double *data, size_t Ni, double *s, double *No, char *flag_anyISNAN, double *W);
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
55 inline int __sumskipnan3w__(double *data, size_t Ni, double *s, double *s2, double *No, char *flag_anyISNAN, double *W);
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
56 inline int __sumskipnan2wr__(double *data, size_t Ni, double *s, double *No, char *flag_anyISNAN, double *W);
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
57 inline int __sumskipnan3wr__(double *data, size_t Ni, double *s, double *s2, double *No, char *flag_anyISNAN, double *W);
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
58 inline int __sumskipnan2we__(double *data, size_t Ni, double *s, double *No, char *flag_anyISNAN, double *W);
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
59 inline int __sumskipnan3we__(double *data, size_t Ni, double *s, double *s2, double *No, char *flag_anyISNAN, double *W);
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
60 inline int __sumskipnan2wer__(double *data, size_t Ni, double *s, double *No, char *flag_anyISNAN, double *W);
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
61 inline int __sumskipnan3wer__(double *data, size_t Ni, double *s, double *s2, double *No, char *flag_anyISNAN, double *W);
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
62
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
63 //#define NO_FLAG
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
64
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
65 #ifdef tmwtypes_h
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
66 #if (MX_API_VER<=0x07020000)
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
67 typedef int mwSize;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
68 #endif
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
69 #endif
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
70
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
71
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
72 void mexFunction(int POutputCount, mxArray* POutput[], int PInputCount, const mxArray *PInputs[])
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
73 {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
74 const mwSize *SZ;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
75 double* LInput;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
76 double* LOutputSum;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
77 double* LOutputCount;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
78 double* LOutputSum2;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
79 long double* LongOutputSum = NULL;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
80 long double* LongOutputCount = NULL;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
81 long double* LongOutputSum2 = NULL;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
82 double x;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
83 double* W = NULL; // weight vector
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
84
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
85 mwSize DIM = 0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
86 mwSize D1, D2, D3; // NN; //
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
87 mwSize ND, ND2; // number of dimensions: input, output
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
88 mwSize ix0, ix1, ix2; // index to input and output
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
89 mwSize j, l; // running indices
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
90 mwSize *SZ2; // size of output
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
91 char flag_isNaN = 0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
92
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
93 // check for proper number of input and output arguments
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
94 if ((PInputCount <= 0) || (PInputCount > 4))
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
95 mexErrMsgTxt("SUMSKIPNAN.MEX requires between 1 and 4 arguments.");
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
96 if (POutputCount > 4)
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
97 mexErrMsgTxt("SUMSKIPNAN.MEX has 1 to 3 output arguments.");
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
98
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
99 // get 1st argument
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
100 if(mxIsDouble(PInputs[0]) && !mxIsComplex(PInputs[0]))
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
101 LInput = mxGetPr(PInputs[0]);
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
102 else
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
103 mexErrMsgTxt("First argument must be REAL/DOUBLE.");
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
104
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
105 // get 2nd argument
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
106 if (PInputCount > 1) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
107 switch (mxGetNumberOfElements(PInputs[1])) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
108 case 0: x = 0.0; // accept empty element
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
109 break;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
110 case 1: x = (mxIsNumeric(PInputs[1]) ? mxGetScalar(PInputs[1]) : -1.0);
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
111 break;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
112 default:x = -1.0; // invalid
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
113 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
114 if ((x < 0) || (x > 65535) || (x != floor(x)))
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
115 mexErrMsgTxt("Error SUMSKIPNAN.MEX: DIM-argument must be a positive integer scalar");
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
116
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
117 DIM = (unsigned)floor(x);
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
118 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
119
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
120 // get size
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
121 ND = mxGetNumberOfDimensions(PInputs[0]);
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
122 // NN = mxGetNumberOfElements(PInputs[0]);
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
123 SZ = mxGetDimensions(PInputs[0]);
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
124
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
125 // if DIM==0 (undefined), look for first dimension with more than 1 element.
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
126 for (j = 0; (DIM < 1) && (j < ND); j++)
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
127 if (SZ[j]>1) DIM = j+1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
128
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
129 if (DIM < 1) DIM=1; // in case DIM is still undefined
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
130
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
131 ND2 = (ND>DIM ? ND : DIM); // number of dimensions of output
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
132
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
133 SZ2 = (mwSize*)mxCalloc(ND2, sizeof(mwSize)); // allocate memory for output size
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
134
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
135 for (j=0; j<ND; j++) // copy size of input;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
136 SZ2[j] = SZ[j];
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
137 for (j=ND; j<ND2; j++) // in case DIM > ND, add extra elements 1
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
138 SZ2[j] = 1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
139
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
140 for (j=0, D1=1; j<DIM-1; D1=D1*SZ2[j++]); // D1 is the number of elements between two elements along dimension DIM
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
141 D2 = SZ2[DIM-1]; // D2 contains the size along dimension DIM
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
142 for (j=DIM, D3=1; j<ND; D3=D3*SZ2[j++]); // D3 is the number of blocks containing D1*D2 elements
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
143
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
144 SZ2[DIM-1] = 1; // size of output is same as size of input but SZ(DIM)=1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
145
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
146 // get weight vector for weighted sumskipnan
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
147 if (PInputCount > 3) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
148 if (!mxGetNumberOfElements(PInputs[3]))
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
149 ; // empty weight vector - no weighting
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
150 else if (mxGetNumberOfElements(PInputs[3])==D2)
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
151 W = mxGetPr(PInputs[3]);
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
152 else
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
153 mexErrMsgTxt("Error SUMSKIPNAN.MEX: length of weight vector does not match size of dimension");
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
154 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
155
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
156 int ACC_LEVEL = 0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
157 {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
158 mxArray *LEVEL = NULL;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
159 int s = mexCallMATLAB(1, &LEVEL, 0, NULL, "flag_accuracy_level");
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
160 if (!s) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
161 ACC_LEVEL = (int) mxGetScalar(LEVEL);
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
162 if ((D1>1) && (ACC_LEVEL>2))
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
163 mexWarnMsgTxt("Warning: Kahan summation not supported with stride > 1 !");
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
164 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
165 mxDestroyArray(LEVEL);
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
166 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
167 // mexPrintf("Accuracy Level=%i\n",ACC_LEVEL);
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
168
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
169 // create outputs
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
170 #define TYP mxDOUBLE_CLASS
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
171
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
172 POutput[0] = mxCreateNumericArray(ND2, SZ2, TYP, mxREAL);
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
173 LOutputSum = mxGetPr(POutput[0]);
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
174 if (D1!=1 && D2>0) LongOutputSum = (long double*) mxCalloc(D1*D3,sizeof(long double));
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
175
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
176 if (POutputCount >= 2) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
177 POutput[1] = mxCreateNumericArray(ND2, SZ2, TYP, mxREAL);
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
178 LOutputCount = mxGetPr(POutput[1]);
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
179 if (D1!=1 && D2>0) LongOutputCount = (long double*) mxCalloc(D1*D3,sizeof(long double));
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
180 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
181 if (POutputCount >= 3) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
182 POutput[2] = mxCreateNumericArray(ND2, SZ2, TYP, mxREAL);
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
183 LOutputSum2 = mxGetPr(POutput[2]);
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
184 if (D1!=1 && D2>0) LongOutputSum2 = (long double*) mxCalloc(D1*D3,sizeof(long double));
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
185 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
186 mxFree(SZ2);
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
187
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
188 if (D1*D2*D3<1) // zero size array
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
189 ; // do nothing
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
190 else if ((D1==1) && (ACC_LEVEL<1)) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
191 // double accuray, naive summation, error = N*2^-52
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
192 switch (POutputCount) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
193 case 1:
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
194 for (l = 0; l<D3; l++) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
195 double count;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
196 __sumskipnan2wr__(LInput+l*D2, D2, LOutputSum+l, &count, &flag_isNaN, W);
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
197 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
198 break;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
199 case 2:
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
200 for (l = 0; l<D3; l++) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
201 __sumskipnan2wr__(LInput+l*D2, D2, LOutputSum+l, LOutputCount+l, &flag_isNaN, W);
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
202 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
203 break;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
204 case 3:
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
205 for (l = 0; l<D3; l++) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
206 __sumskipnan3wr__(LInput+l*D2, D2, LOutputSum+l, LOutputSum2+l, LOutputCount+l, &flag_isNaN, W);
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
207 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
208 break;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
209 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
210 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
211 else if ((D1==1) && (ACC_LEVEL==1)) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
212 // extended accuray, naive summation, error = N*2^-64
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
213 switch (POutputCount) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
214 case 1:
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
215 for (l = 0; l<D3; l++) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
216 double count;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
217 __sumskipnan2w__(LInput+l*D2, D2, LOutputSum+l, &count, &flag_isNaN, W);
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
218 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
219 break;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
220 case 2:
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
221 for (l = 0; l<D3; l++) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
222 __sumskipnan2w__(LInput+l*D2, D2, LOutputSum+l, LOutputCount+l, &flag_isNaN, W);
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
223 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
224 break;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
225 case 3:
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
226 for (l = 0; l<D3; l++) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
227 __sumskipnan3w__(LInput+l*D2, D2, LOutputSum+l, LOutputSum2+l, LOutputCount+l, &flag_isNaN, W);
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
228 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
229 break;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
230 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
231 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
232 else if ((D1==1) && (ACC_LEVEL==3)) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
233 // ACC_LEVEL==3: extended accuracy and Kahan Summation, error = 2^-64
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
234 switch (POutputCount) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
235 case 1:
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
236 for (l = 0; l<D3; l++) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
237 double count;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
238 __sumskipnan2we__(LInput+l*D2, D2, LOutputSum+l, &count, &flag_isNaN, W);
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
239 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
240 break;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
241 case 2:
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
242 for (l = 0; l<D3; l++) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
243 __sumskipnan2we__(LInput+l*D2, D2, LOutputSum+l, LOutputCount+l, &flag_isNaN, W);
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
244 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
245 break;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
246 case 3:
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
247 for (l = 0; l<D3; l++) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
248 __sumskipnan3we__(LInput+l*D2, D2, LOutputSum+l, LOutputSum2+l, LOutputCount+l, &flag_isNaN, W);
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
249 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
250 break;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
251 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
252 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
253 else if ((D1==1) && (ACC_LEVEL==2)) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
254 // ACC_LEVEL==2: double accuracy and Kahan Summation, error = 2^-52
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
255 switch (POutputCount) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
256 case 1:
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
257 for (l = 0; l<D3; l++) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
258 double count;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
259 __sumskipnan2wer__(LInput+l*D2, D2, LOutputSum+l, &count, &flag_isNaN, W);
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
260 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
261 break;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
262 case 2:
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
263 for (l = 0; l<D3; l++) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
264 __sumskipnan2wer__(LInput+l*D2, D2, LOutputSum+l, LOutputCount+l, &flag_isNaN, W);
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
265 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
266 break;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
267 case 3:
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
268 for (l = 0; l<D3; l++) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
269 __sumskipnan3wer__(LInput+l*D2, D2, LOutputSum+l, LOutputSum2+l, LOutputCount+l, &flag_isNaN, W);
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
270 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
271 break;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
272 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
273 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
274 else if (POutputCount <= 1) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
275 // OUTER LOOP: along dimensions > DIM
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
276 for (l = 0; l<D3; l++) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
277 ix0 = l*D1; // index for output
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
278 ix1 = ix0*D2; // index for input
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
279 for (j=0; j<D2; j++) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
280 // minimize cache misses
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
281 ix2 = ix0; // index for output
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
282 // Inner LOOP: along dimensions < DIM
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
283 if (W) do {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
284 long double x = *LInput;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
285 if (!isnan(x)) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
286 LongOutputSum[ix2] += W[j]*x;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
287 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
288 #ifndef NO_FLAG
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
289 else
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
290 flag_isNaN = 1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
291 #endif
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
292 LInput++;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
293 ix2++;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
294 } while (ix2 != (l+1)*D1);
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
295 else do {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
296 long double x = *LInput;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
297 if (!isnan(x)) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
298 LongOutputSum[ix2] += x;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
299 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
300 #ifndef NO_FLAG
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
301 else
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
302 flag_isNaN = 1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
303 #endif
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
304 LInput++;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
305 ix2++;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
306 } while (ix2 != (l+1)*D1);
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
307 } // end for (j=
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
308
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
309 /* copy to output */
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
310 for (j=0; j<D1; j++) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
311 LOutputSum[ix0+j] = LongOutputSum[ix0+j];
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
312 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
313 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
314 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
315
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
316 else if (POutputCount == 2) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
317 // OUTER LOOP: along dimensions > DIM
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
318 for (l = 0; l<D3; l++) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
319 ix0 = l*D1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
320 ix1 = ix0*D2; // index for input
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
321 for (j=0; j<D2; j++) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
322 // minimize cache misses
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
323 ix2 = ix0; // index for output
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
324 // Inner LOOP: along dimensions < DIM
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
325 if (W) do {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
326 long double x = *LInput;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
327 if (!isnan(x)) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
328 LongOutputCount[ix2] += W[j];
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
329 LongOutputSum[ix2] += W[j]*x;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
330 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
331 #ifndef NO_FLAG
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
332 else
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
333 flag_isNaN = 1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
334 #endif
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
335 LInput++;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
336 ix2++;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
337 } while (ix2 != (l+1)*D1);
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
338 else do {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
339 long double x = *LInput;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
340 if (!isnan(x)) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
341 LongOutputCount[ix2] += 1.0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
342 LongOutputSum[ix2] += x;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
343 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
344 #ifndef NO_FLAG
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
345 else
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
346 flag_isNaN = 1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
347 #endif
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
348 LInput++;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
349 ix2++;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
350 } while (ix2 != (l+1)*D1);
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
351 } // end for (j=
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
352
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
353 /* copy to output */
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
354 for (j=0; j<D1; j++) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
355 LOutputSum[ix0+j] = LongOutputSum[ix0+j];
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
356 LOutputCount[ix0+j] = LongOutputCount[ix0+j];
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
357 } // end else
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
358 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
359 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
360
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
361 else if (POutputCount == 3) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
362 // OUTER LOOP: along dimensions > DIM
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
363 for (l = 0; l<D3; l++) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
364 ix0 = l*D1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
365 ix1 = ix0*D2; // index for input
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
366 for (j=0; j<D2; j++) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
367 // minimize cache misses
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
368 ix2 = ix0; // index for output
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
369 // Inner LOOP: along dimensions < DIM
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
370 if (W) do {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
371 long double x = *LInput;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
372 if (!isnan(x)) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
373 LongOutputCount[ix2] += W[j];
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
374 double t = W[j]*x;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
375 LongOutputSum[ix2] += t;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
376 LongOutputSum2[ix2] += x*t;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
377 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
378 #ifndef NO_FLAG
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
379 else
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
380 flag_isNaN = 1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
381 #endif
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
382 LInput++;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
383 ix2++;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
384 } while (ix2 != (l+1)*D1);
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
385 else do {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
386 long double x = *LInput;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
387 if (!isnan(x)) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
388 LongOutputCount[ix2] += 1.0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
389 LongOutputSum[ix2] += x;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
390 LongOutputSum2[ix2] += x*x;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
391 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
392 #ifndef NO_FLAG
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
393 else
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
394 flag_isNaN = 1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
395 #endif
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
396 LInput++;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
397 ix2++;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
398 } while (ix2 != (l+1)*D1);
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
399 } // end for (j=
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
400
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
401 /* copy to output */
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
402 for (j=0; j<D1; j++) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
403 LOutputSum[ix0+j] = LongOutputSum[ix0+j];
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
404 LOutputCount[ix0+j] = LongOutputCount[ix0+j];
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
405 LOutputSum2[ix0+j] = LongOutputSum2[ix0+j];
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
406 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
407 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
408 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
409 if (LongOutputSum) mxFree(LongOutputSum);
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
410 if (LongOutputCount) mxFree(LongOutputCount);
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
411 if (LongOutputSum2) mxFree(LongOutputSum2);
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
412
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
413 #ifndef NO_FLAG
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
414 //mexPrintf("Third argument must be not empty - otherwise status whether a NaN occured or not cannot be returned.");
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
415 /* this is a hack, the third input argument is used to return whether a NaN occured or not.
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
416 this requires that the input argument is a non-empty variable
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
417 */
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
418 if (flag_isNaN && (PInputCount > 2) && mxGetNumberOfElements(PInputs[2])) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
419 // set FLAG_NANS_OCCURED
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
420 switch (mxGetClassID(PInputs[2])) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
421 case mxLOGICAL_CLASS:
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
422 case mxCHAR_CLASS:
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
423 case mxINT8_CLASS:
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
424 case mxUINT8_CLASS:
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
425 *(uint8_t*)mxGetData(PInputs[2]) = 1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
426 break;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
427 case mxDOUBLE_CLASS:
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
428 *(double*)mxGetData(PInputs[2]) = 1.0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
429 break;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
430 case mxSINGLE_CLASS:
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
431 *(float*)mxGetData(PInputs[2]) = 1.0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
432 break;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
433 case mxINT16_CLASS:
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
434 case mxUINT16_CLASS:
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
435 *(uint16_t*)mxGetData(PInputs[2]) = 1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
436 break;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
437 case mxINT32_CLASS:
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
438 case mxUINT32_CLASS:
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
439 *(uint32_t*)mxGetData(PInputs[2])= 1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
440 break;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
441 case mxINT64_CLASS:
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
442 case mxUINT64_CLASS:
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
443 *(uint64_t*)mxGetData(PInputs[2]) = 1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
444 break;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
445 case mxFUNCTION_CLASS:
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
446 case mxUNKNOWN_CLASS:
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
447 case mxCELL_CLASS:
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
448 case mxSTRUCT_CLASS:
7301
485d1594d155 addresses undesired side-effect off in-place sorting of data
schloegl
parents: 6549
diff changeset
449 default:
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
450 mexPrintf("Type of 3rd input argument not supported.");
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
451 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
452 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
453 #endif
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
454 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
455
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
456
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
457 #define stride 1
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
458 inline int __sumskipnan2w__(double *data, size_t Ni, double *s, double *No, char *flag_anyISNAN, double *W)
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
459 {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
460 long double sum=0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
461 char flag=0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
462 // LOOP along dimension DIM
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
463
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
464 double *end = data + stride*Ni;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
465 if (W) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
466 // with weight vector
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
467 long double count = 0.0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
468 do {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
469 long double x = *data;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
470 if (!isnan(x))
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
471 {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
472 count += *W;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
473 sum += *W*x;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
474 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
475 #ifndef NO_FLAG
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
476 else
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
477 flag = 1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
478 #endif
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
479
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
480 data++; // stride=1
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
481 W++;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
482 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
483 while (data < end);
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
484 *No = count;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
485 } else {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
486 // w/o weight vector
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
487 size_t countI = 0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
488 do {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
489 long double x = *data;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
490 if (!isnan(x))
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
491 {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
492 countI++;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
493 sum += x;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
494 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
495 #ifndef NO_FLAG
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
496 else
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
497 flag = 1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
498 #endif
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
499 data++; // stride=1
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
500 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
501 while (data < end);
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
502 *No = (double)countI;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
503 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
504
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
505 #ifndef NO_FLAG
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
506 if (flag && (flag_anyISNAN != NULL)) *flag_anyISNAN = 1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
507 #endif
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
508 *s = sum;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
509
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
510 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
511
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
512
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
513 inline int __sumskipnan3w__(double *data, size_t Ni, double *s, double *s2, double *No, char *flag_anyISNAN, double *W)
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
514 {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
515 long double sum=0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
516 long double msq=0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
517 char flag=0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
518 // LOOP along dimension DIM
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
519
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
520 double *end = data + stride*Ni;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
521 if (W) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
522 // with weight vector
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
523 long double count = 0.0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
524 do {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
525 long double x = *data;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
526 if (!isnan(x)) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
527 count += *W;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
528 long double t = *W*x;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
529 sum += t;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
530 msq += x*t;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
531 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
532 #ifndef NO_FLAG
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
533 else
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
534 flag = 1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
535 #endif
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
536 data++; // stride=1
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
537 W++;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
538 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
539 while (data < end);
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
540 *No = count;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
541 } else {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
542 // w/o weight vector
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
543 size_t countI = 0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
544 do {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
545 long double x = *data;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
546 if (!isnan(x)) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
547 countI++;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
548 sum += x;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
549 msq += x*x;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
550 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
551 #ifndef NO_FLAG
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
552 else
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
553 flag = 1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
554 #endif
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
555 data++; // stride=1
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
556 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
557 while (data < end);
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
558 *No = (double)countI;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
559 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
560
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
561 #ifndef NO_FLAG
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
562 if (flag && (flag_anyISNAN != NULL)) *flag_anyISNAN = 1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
563 #endif
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
564 *s = sum;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
565 *s2 = msq;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
566 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
567
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
568 inline int __sumskipnan2wr__(double *data, size_t Ni, double *s, double *No, char *flag_anyISNAN, double *W)
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
569 {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
570 double sum=0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
571 char flag=0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
572 // LOOP along dimension DIM
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
573
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
574 double *end = data + stride*Ni;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
575 if (W) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
576 // with weight vector
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
577 double count = 0.0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
578 do {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
579 double x = *data;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
580 if (!isnan(x))
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
581 {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
582 count += *W;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
583 sum += *W*x;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
584 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
585 #ifndef NO_FLAG
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
586 else
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
587 flag = 1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
588 #endif
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
589
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
590 data++; // stride=1
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
591 W++;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
592 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
593 while (data < end);
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
594 *No = count;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
595 } else {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
596 // w/o weight vector
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
597 size_t countI = 0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
598 do {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
599 double x = *data;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
600 if (!isnan(x))
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
601 {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
602 countI++;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
603 sum += x;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
604 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
605 #ifndef NO_FLAG
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
606 else
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
607 flag = 1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
608 #endif
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
609 data++; // stride=1
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
610 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
611 while (data < end);
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
612 *No = (double)countI;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
613 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
614
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
615 #ifndef NO_FLAG
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
616 if (flag && (flag_anyISNAN != NULL)) *flag_anyISNAN = 1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
617 #endif
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
618 *s = sum;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
619
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
620 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
621
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
622
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
623 inline int __sumskipnan3wr__(double *data, size_t Ni, double *s, double *s2, double *No, char *flag_anyISNAN, double *W)
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
624 {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
625 double sum=0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
626 double msq=0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
627 char flag=0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
628 // LOOP along dimension DIM
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
629
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
630 double *end = data + stride*Ni;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
631 if (W) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
632 // with weight vector
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
633 double count = 0.0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
634 do {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
635 double x = *data;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
636 if (!isnan(x)) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
637 count += *W;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
638 double t = *W*x;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
639 sum += t;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
640 msq += x*t;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
641 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
642 #ifndef NO_FLAG
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
643 else
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
644 flag = 1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
645 #endif
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
646 data++; // stride=1
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
647 W++;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
648 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
649 while (data < end);
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
650 *No = count;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
651 } else {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
652 // w/o weight vector
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
653 size_t countI = 0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
654 do {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
655 double x = *data;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
656 if (!isnan(x)) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
657 countI++;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
658 sum += x;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
659 msq += x*x;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
660 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
661 #ifndef NO_FLAG
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
662 else
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
663 flag = 1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
664 #endif
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
665 data++; // stride=1
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
666 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
667 while (data < end);
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
668 *No = (double)countI;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
669 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
670
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
671 #ifndef NO_FLAG
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
672 if (flag && (flag_anyISNAN != NULL)) *flag_anyISNAN = 1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
673 #endif
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
674 *s = sum;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
675 *s2 = msq;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
676 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
677
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
678
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
679
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
680 /***************************************
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
681 using Kahan's summation formula [1]
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
682 this gives more accurate results while the computational effort within the loop is about 4x as high
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
683 First tests show a penalty of about 40% in terms of computational time.
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
684
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
685 [1] David Goldberg,
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
686 What Every Computer Scientist Should Know About Floating-Point Arithmetic
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
687 ACM Computing Surveys, Vol 23, No 1, March 1991.
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
688 ****************************************/
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
689
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
690 inline int __sumskipnan2we__(double *data, size_t Ni, double *s, double *No, char *flag_anyISNAN, double *W)
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
691 {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
692 long double sum=0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
693 char flag=0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
694 // LOOP along dimension DIM
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
695
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
696 double *end = data + stride*Ni;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
697 if (W) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
698 // with weight vector
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
699 long double count = 0.0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
700 long double rc=0.0, rn=0.0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
701 do {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
702 long double x = *data;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
703 long double t,y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
704 if (!isnan(x))
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
705 {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
706 //count += *W; [1]
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
707 y = *W-rn;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
708 t = count+y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
709 rn= (t-count)-y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
710 count= t;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
711
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
712 //sum += *W*x; [1]
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
713 y = *W*x-rc;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
714 t = sum+y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
715 rc= (t-sum)-y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
716 sum= t;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
717 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
718 #ifndef NO_FLAG
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
719 else
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
720 flag = 1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
721 #endif
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
722
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
723 data++; // stride=1
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
724 W++;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
725 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
726 while (data < end);
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
727 *No = count;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
728 } else {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
729 // w/o weight vector
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
730 size_t countI = 0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
731 long double rc=0.0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
732 do {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
733 long double x = *data;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
734 long double t,y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
735 if (!isnan(x))
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
736 {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
737 countI++;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
738 // sum += x; [1]
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
739 y = x-rc;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
740 t = sum+y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
741 rc= (t-sum)-y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
742 sum= t;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
743 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
744 #ifndef NO_FLAG
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
745 else
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
746 flag = 1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
747 #endif
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
748 data++; // stride=1
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
749 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
750 while (data < end);
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
751 *No = (double)countI;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
752 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
753
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
754 #ifndef NO_FLAG
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
755 if (flag && (flag_anyISNAN != NULL)) *flag_anyISNAN = 1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
756 #endif
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
757 *s = sum;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
758
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
759 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
760
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
761
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
762 inline int __sumskipnan3we__(double *data, size_t Ni, double *s, double *s2, double *No, char *flag_anyISNAN, double *W)
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
763 {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
764 long double sum=0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
765 long double msq=0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
766 char flag=0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
767 // LOOP along dimension DIM
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
768
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
769 double *end = data + stride*Ni;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
770 if (W) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
771 // with weight vector
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
772 long double count = 0.0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
773 long double rc=0.0, rn=0.0, rq=0.0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
774 do {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
775 long double x = *data;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
776 long double t,y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
777 if (!isnan(x)) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
778 //count += *W; [1]
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
779 y = *W-rn;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
780 t = count+y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
781 rn= (t-count)-y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
782 count= t;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
783
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
784 long double w = *W*x;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
785 //sum += *W*x; [1]
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
786 y = *W*x-rc;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
787 t = sum+y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
788 rc= (t-sum)-y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
789 sum= t;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
790
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
791 // msq += x*w;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
792 y = w*x-rq;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
793 t = msq+y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
794 rq= (t-msq)-y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
795 msq= t;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
796 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
797 #ifndef NO_FLAG
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
798 else
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
799 flag = 1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
800 #endif
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
801 data++; // stride=1
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
802 W++;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
803 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
804 while (data < end);
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
805 *No = count;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
806 } else {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
807 // w/o weight vector
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
808 size_t countI = 0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
809 long double rc=0.0, rq=0.0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
810 do {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
811 long double x = *data;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
812 long double t,y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
813 if (!isnan(x)) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
814 countI++;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
815 //sum += x; [1]
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
816 y = x-rc;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
817 t = sum+y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
818 rc= (t-sum)-y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
819 sum= t;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
820
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
821 // msq += x*x;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
822 y = x*x-rq;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
823 t = msq+y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
824 rq= (t-msq)-y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
825 msq= t;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
826 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
827 #ifndef NO_FLAG
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
828 else
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
829 flag = 1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
830 #endif
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
831 data++; // stride=1
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
832 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
833 while (data < end);
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
834 *No = (double)countI;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
835 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
836
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
837 #ifndef NO_FLAG
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
838 if (flag && (flag_anyISNAN != NULL)) *flag_anyISNAN = 1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
839 #endif
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
840 *s = sum;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
841 *s2 = msq;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
842 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
843
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
844 inline int __sumskipnan2wer__(double *data, size_t Ni, double *s, double *No, char *flag_anyISNAN, double *W)
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
845 {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
846 double sum=0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
847 char flag=0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
848 // LOOP along dimension DIM
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
849
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
850 double *end = data + stride*Ni;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
851 if (W) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
852 // with weight vector
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
853 double count = 0.0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
854 double rc=0.0, rn=0.0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
855 do {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
856 double x = *data;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
857 double t,y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
858 if (!isnan(x))
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
859 {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
860 //count += *W; [1]
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
861 y = *W-rn;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
862 t = count+y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
863 rn= (t-count)-y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
864 count= t;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
865
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
866 //sum += *W*x; [1]
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
867 y = *W*x-rc;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
868 t = sum+y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
869 rc= (t-sum)-y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
870 sum= t;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
871 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
872 #ifndef NO_FLAG
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
873 else
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
874 flag = 1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
875 #endif
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
876
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
877 data++; // stride=1
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
878 W++;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
879 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
880 while (data < end);
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
881 *No = count;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
882 } else {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
883 // w/o weight vector
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
884 size_t countI = 0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
885 double rc=0.0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
886 do {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
887 double x = *data;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
888 double t,y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
889 if (!isnan(x))
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
890 {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
891 countI++;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
892 // sum += x; [1]
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
893 y = x-rc;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
894 t = sum+y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
895 rc= (t-sum)-y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
896 sum= t;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
897 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
898 #ifndef NO_FLAG
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
899 else
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
900 flag = 1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
901 #endif
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
902 data++; // stride=1
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
903 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
904 while (data < end);
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
905 *No = (double)countI;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
906 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
907
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
908 #ifndef NO_FLAG
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
909 if (flag && (flag_anyISNAN != NULL)) *flag_anyISNAN = 1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
910 #endif
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
911 *s = sum;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
912
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
913 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
914
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
915
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
916 inline int __sumskipnan3wer__(double *data, size_t Ni, double *s, double *s2, double *No, char *flag_anyISNAN, double *W)
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
917 {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
918 double sum=0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
919 double msq=0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
920 char flag=0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
921 // LOOP along dimension DIM
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
922
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
923 double *end = data + stride*Ni;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
924 if (W) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
925 // with weight vector
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
926 double count = 0.0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
927 double rc=0.0, rn=0.0, rq=0.0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
928 do {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
929 double x = *data;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
930 double t,y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
931 if (!isnan(x)) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
932 //count += *W; [1]
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
933 y = *W-rn;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
934 t = count+y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
935 rn= (t-count)-y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
936 count= t;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
937
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
938 double w = *W*x;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
939 //sum += *W*x; [1]
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
940 y = *W*x-rc;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
941 t = sum+y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
942 rc= (t-sum)-y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
943 sum= t;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
944
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
945 // msq += x*w;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
946 y = w*x-rq;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
947 t = msq+y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
948 rq= (t-msq)-y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
949 msq= t;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
950 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
951 #ifndef NO_FLAG
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
952 else
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
953 flag = 1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
954 #endif
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
955 data++; // stride=1
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
956 W++;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
957 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
958 while (data < end);
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
959 *No = count;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
960 } else {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
961 // w/o weight vector
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
962 size_t countI = 0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
963 double rc=0.0, rq=0.0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
964 do {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
965 double x = *data;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
966 double t,y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
967 if (!isnan(x)) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
968 countI++;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
969 //sum += x; [1]
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
970 y = x-rc;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
971 t = sum+y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
972 rc= (t-sum)-y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
973 sum= t;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
974
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
975 // msq += x*x;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
976 y = x*x-rq;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
977 t = msq+y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
978 rq= (t-msq)-y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
979 msq= t;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
980 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
981 #ifndef NO_FLAG
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
982 else
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
983 flag = 1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
984 #endif
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
985 data++; // stride=1
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
986 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
987 while (data < end);
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
988 *No = (double)countI;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
989 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
990
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
991 #ifndef NO_FLAG
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
992 if (flag && (flag_anyISNAN != NULL)) *flag_anyISNAN = 1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
993 #endif
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
994 *s = sum;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
995 *s2 = msq;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
996 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
997