Mercurial > forge

annotate extra/NaN/src/sumskipnan_mex.cpp @ 6549:41e9854fe26d octave-forge

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