Mercurial > forge

comparison extra/NaN/src/covm_mex.cpp @ 12691:6d6285a2a633 octave-forge

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
use macro ISNAN() instead of C++'s isnan() - because it supports all floating point formats not just double
author schloegl
date Sat, 12 Sep 2015 14:16:39 +0000
parents f26b1170ea90
children 79e7259c6ff1
comparison
equal deleted inserted replaced
12690:16bc2657b1f1 12691:6d6285a2a633
27 // - flag: if not empty, it is set to 1 if some NaN was observed 27 // - flag: if not empty, it is set to 1 if some NaN was observed
28 // - W: weight vector to compute weighted correlation 28 // - W: weight vector to compute weighted correlation
29 // 29 //
30 // Output: 30 // Output:
31 // - CC = X' * sparse(diag(W)) * Y while NaN's are skipped 31 // - CC = X' * sparse(diag(W)) * Y while NaN's are skipped
32 // - NN = real(~isnan(X)')*sparse(diag(W))*real(~isnan(Y)) count of valid (non-NaN) elements 32 // - NN = real(~ISNAN(X)')*sparse(diag(W))*real(~ISNAN(Y)) count of valid (non-NaN) elements
33 // computed more efficiently 33 // computed more efficiently
34 // 34 //
35 // $Id$ 35 // $Id$
36 // Copyright (C) 2009,2010,2011 Alois Schloegl <alois.schloegl@gmail.com> 36 // Copyright (C) 2009,2010,2011 Alois Schloegl <alois.schloegl@gmail.com>
37 // This function is part of the NaN-toolbox 37 // This function is part of the NaN-toolbox
46 #include <math.h> 46 #include <math.h>
47 #include "mex.h" 47 #include "mex.h"
48 48
49 /*#define NO_FLAG*/ 49 /*#define NO_FLAG*/
50 50
51 /*
52 math.h has isnan() defined for all sizes of floating point numbers,
53 but c++ assumes isnan(double), causing possible conversions for float and long double
54 */
55 #define ISNAN(a) (a!=a)
51 56
52 void mexFunction(int POutputCount, mxArray* POutput[], int PInputCount, const mxArray *PInputs[]) 57 void mexFunction(int POutputCount, mxArray* POutput[], int PInputCount, const mxArray *PInputs[])
53 { 58 {
54 double *X0=NULL, *Y0=NULL, *W=NULL; 59 double *X0=NULL, *Y0=NULL, *W=NULL;
55 double *CC; 60 double *CC;
161 double w; 166 double w;
162 if (W) { 167 if (W) {
163 w = W[k]; 168 w = W[k];
164 for (i=0; i<cX; i++) { 169 for (i=0; i<cX; i++) {
165 double x = X0[k+i*rX]; 170 double x = X0[k+i*rX];
166 if (isnan(x)) { 171 if (ISNAN(x)) {
167 #ifndef NO_FLAG 172 #ifndef NO_FLAG
168 flag_isNaN = 1; 173 flag_isNaN = 1;
169 #endif 174 #endif
170 continue; 175 continue;
171 } 176 }
172 for (j=0; j<cY; j++) { 177 for (j=0; j<cY; j++) {
173 double y = Y0[k+j*rY]; 178 double y = Y0[k+j*rY];
174 if (isnan(y)) { 179 if (ISNAN(y)) {
175 #ifndef NO_FLAG 180 #ifndef NO_FLAG
176 flag_isNaN = 1; 181 flag_isNaN = 1;
177 #endif 182 #endif
178 continue; 183 continue;
179 } 184 }
183 } 188 }
184 } 189 }
185 } 190 }
186 else for (i=0; i<cX; i++) { 191 else for (i=0; i<cX; i++) {
187 double x = X0[k+i*rX]; 192 double x = X0[k+i*rX];
188 if (isnan(x)) { 193 if (ISNAN(x)) {
189 #ifndef NO_FLAG 194 #ifndef NO_FLAG
190 flag_isNaN = 1; 195 flag_isNaN = 1;
191 #endif 196 #endif
192 continue; 197 continue;
193 } 198 }
194 for (j=0; j<cY; j++) { 199 for (j=0; j<cY; j++) {
195 double y = Y0[k+j*rY]; 200 double y = Y0[k+j*rY];
196 if (isnan(y)) { 201 if (ISNAN(y)) {
197 #ifndef NO_FLAG 202 #ifndef NO_FLAG
198 flag_isNaN = 1; 203 flag_isNaN = 1;
199 #endif 204 #endif
200 continue; 205 continue;
201 } 206 }
228 double cc = 0.0; 233 double cc = 0.0;
229 double nw = 0.0; 234 double nw = 0.0;
230 size_t k; 235 size_t k;
231 for (k=0; k<rX; k++) { 236 for (k=0; k<rX; k++) {
232 double z = X[k]*Y[k]; 237 double z = X[k]*Y[k];
233 if (isnan(z)) { 238 if (ISNAN(z)) {
234 #ifndef NO_FLAG 239 #ifndef NO_FLAG
235 flag_isNaN = 1; 240 flag_isNaN = 1;
236 #endif 241 #endif
237 continue; 242 continue;
238 } 243 }
252 double cc = 0.0; 257 double cc = 0.0;
253 size_t nn = 0; 258 size_t nn = 0;
254 size_t k; 259 size_t k;
255 for (k=0; k<rX; k++) { 260 for (k=0; k<rX; k++) {
256 double z = X[k]*Y[k]; 261 double z = X[k]*Y[k];
257 if (isnan(z)) { 262 if (ISNAN(z)) {
258 #ifndef NO_FLAG 263 #ifndef NO_FLAG
259 flag_isNaN = 1; 264 flag_isNaN = 1;
260 #endif 265 #endif
261 continue; 266 continue;
262 } 267 }
281 double cc = 0.0; 286 double cc = 0.0;
282 double nw = 0.0; 287 double nw = 0.0;
283 size_t k; 288 size_t k;
284 for (k=0; k<rX; k++) { 289 for (k=0; k<rX; k++) {
285 double z = X[k]*Y[k]; 290 double z = X[k]*Y[k];
286 if (isnan(z)) { 291 if (ISNAN(z)) {
287 #ifndef NO_FLAG 292 #ifndef NO_FLAG
288 flag_isNaN = 1; 293 flag_isNaN = 1;
289 #endif 294 #endif
290 continue; 295 continue;
291 } 296 }
312 double cc = 0.0; 317 double cc = 0.0;
313 size_t nn = 0; 318 size_t nn = 0;
314 size_t k; 319 size_t k;
315 for (k=0; k<rX; k++) { 320 for (k=0; k<rX; k++) {
316 double z = X[k]*Y[k]; 321 double z = X[k]*Y[k];
317 if (isnan(z)) { 322 if (ISNAN(z)) {
318 #ifndef NO_FLAG 323 #ifndef NO_FLAG
319 flag_isNaN = 1; 324 flag_isNaN = 1;
320 #endif 325 #endif
321 continue; 326 continue;
322 } 327 }
351 long double cc=0.0; 356 long double cc=0.0;
352 long double nn=0.0; 357 long double nn=0.0;
353 size_t k; 358 size_t k;
354 for (k=0; k<rX; k++) { 359 for (k=0; k<rX; k++) {
355 long double z = ((long double)X[k])*Y[k]; 360 long double z = ((long double)X[k])*Y[k];
356 if (isnan(z)) { 361 if (ISNAN(z)) {
357 #ifndef NO_FLAG 362 #ifndef NO_FLAG
358 flag_isNaN = 1; 363 flag_isNaN = 1;
359 #endif 364 #endif
360 continue; 365 continue;
361 } 366 }
375 long double cc=0.0; 380 long double cc=0.0;
376 size_t nn=0; 381 size_t nn=0;
377 size_t k; 382 size_t k;
378 for (k=0; k<rX; k++) { 383 for (k=0; k<rX; k++) {
379 long double z = ((long double)X[k])*Y[k]; 384 long double z = ((long double)X[k])*Y[k];
380 if (isnan(z)) { 385 if (ISNAN(z)) {
381 #ifndef NO_FLAG 386 #ifndef NO_FLAG
382 flag_isNaN = 1; 387 flag_isNaN = 1;
383 #endif 388 #endif
384 continue; 389 continue;
385 } 390 }
404 long double cc=0.0; 409 long double cc=0.0;
405 long double nn=0.0; 410 long double nn=0.0;
406 size_t k; 411 size_t k;
407 for (k=0; k<rX; k++) { 412 for (k=0; k<rX; k++) {
408 long double z = ((long double)X[k])*Y[k]; 413 long double z = ((long double)X[k])*Y[k];
409 if (isnan(z)) { 414 if (ISNAN(z)) {
410 #ifndef NO_FLAG 415 #ifndef NO_FLAG
411 flag_isNaN = 1; 416 flag_isNaN = 1;
412 #endif 417 #endif
413 continue; 418 continue;
414 } 419 }
435 long double cc=0.0; 440 long double cc=0.0;
436 size_t nn=0; 441 size_t nn=0;
437 size_t k; 442 size_t k;
438 for (k=0; k<rX; k++) { 443 for (k=0; k<rX; k++) {
439 long double z = ((long double)X[k])*Y[k]; 444 long double z = ((long double)X[k])*Y[k];
440 if (isnan(z)) { 445 if (ISNAN(z)) {
441 #ifndef NO_FLAG 446 #ifndef NO_FLAG
442 flag_isNaN = 1; 447 flag_isNaN = 1;
443 #endif 448 #endif
444 continue; 449 continue;
445 } 450 }
480 long double rn=0.0; 485 long double rn=0.0;
481 size_t k; 486 size_t k;
482 for (k=0; k<rX; k++) { 487 for (k=0; k<rX; k++) {
483 long double t,y; 488 long double t,y;
484 long double z = ((long double)X[k])*Y[k]; 489 long double z = ((long double)X[k])*Y[k];
485 if (isnan(z)) { 490 if (ISNAN(z)) {
486 #ifndef NO_FLAG 491 #ifndef NO_FLAG
487 flag_isNaN = 1; 492 flag_isNaN = 1;
488 #endif 493 #endif
489 continue; 494 continue;
490 } 495 }
515 size_t nn=0; 520 size_t nn=0;
516 size_t k; 521 size_t k;
517 for (k=0; k<rX; k++) { 522 for (k=0; k<rX; k++) {
518 long double t,y; 523 long double t,y;
519 long double z = ((long double)X[k])*Y[k]; 524 long double z = ((long double)X[k])*Y[k];
520 if (isnan(z)) { 525 if (ISNAN(z)) {
521 #ifndef NO_FLAG 526 #ifndef NO_FLAG
522 flag_isNaN = 1; 527 flag_isNaN = 1;
523 #endif 528 #endif
524 continue; 529 continue;
525 } 530 }
552 long double rn=0.0; 557 long double rn=0.0;
553 size_t k; 558 size_t k;
554 for (k=0; k<rX; k++) { 559 for (k=0; k<rX; k++) {
555 long double t,y; 560 long double t,y;
556 long double z = ((long double)X[k])*Y[k]; 561 long double z = ((long double)X[k])*Y[k];
557 if (isnan(z)) { 562 if (ISNAN(z)) {
558 #ifndef NO_FLAG 563 #ifndef NO_FLAG
559 flag_isNaN = 1; 564 flag_isNaN = 1;
560 #endif 565 #endif
561 continue; 566 continue;
562 } 567 }
594 size_t nn=0; 599 size_t nn=0;
595 size_t k; 600 size_t k;
596 for (k=0; k<rX; k++) { 601 for (k=0; k<rX; k++) {
597 long double t,y; 602 long double t,y;
598 long double z = ((long double)X[k])*Y[k]; 603 long double z = ((long double)X[k])*Y[k];
599 if (isnan(z)) { 604 if (ISNAN(z)) {
600 #ifndef NO_FLAG 605 #ifndef NO_FLAG
601 flag_isNaN = 1; 606 flag_isNaN = 1;
602 #endif 607 #endif
603 continue; 608 continue;
604 } 609 }
643 double rn=0.0; 648 double rn=0.0;
644 size_t k; 649 size_t k;
645 for (k=0; k<rX; k++) { 650 for (k=0; k<rX; k++) {
646 double t,y; 651 double t,y;
647 double z = X[k]*Y[k]; 652 double z = X[k]*Y[k];
648 if (isnan(z)) { 653 if (ISNAN(z)) {
649 #ifndef NO_FLAG 654 #ifndef NO_FLAG
650 flag_isNaN = 1; 655 flag_isNaN = 1;
651 #endif 656 #endif
652 continue; 657 continue;
653 } 658 }
678 size_t nn=0; 683 size_t nn=0;
679 size_t k; 684 size_t k;
680 for (k=0; k<rX; k++) { 685 for (k=0; k<rX; k++) {
681 double t,y; 686 double t,y;
682 double z = X[k]*Y[k]; 687 double z = X[k]*Y[k];
683 if (isnan(z)) { 688 if (ISNAN(z)) {
684 #ifndef NO_FLAG 689 #ifndef NO_FLAG
685 flag_isNaN = 1; 690 flag_isNaN = 1;
686 #endif 691 #endif
687 continue; 692 continue;
688 } 693 }
715 double rn=0.0; 720 double rn=0.0;
716 size_t k; 721 size_t k;
717 for (k=0; k<rX; k++) { 722 for (k=0; k<rX; k++) {
718 double t,y; 723 double t,y;
719 double z = X[k]*Y[k]; 724 double z = X[k]*Y[k];
720 if (isnan(z)) { 725 if (ISNAN(z)) {
721 #ifndef NO_FLAG 726 #ifndef NO_FLAG
722 flag_isNaN = 1; 727 flag_isNaN = 1;
723 #endif 728 #endif
724 continue; 729 continue;
725 } 730 }
757 size_t nn=0; 762 size_t nn=0;
758 size_t k; 763 size_t k;
759 for (k=0; k<rX; k++) { 764 for (k=0; k<rX; k++) {
760 double t,y; 765 double t,y;
761 double z = X[k]*Y[k]; 766 double z = X[k]*Y[k];
762 if (isnan(z)) { 767 if (ISNAN(z)) {
763 #ifndef NO_FLAG 768 #ifndef NO_FLAG
764 flag_isNaN = 1; 769 flag_isNaN = 1;
765 #endif 770 #endif
766 continue; 771 continue;
767 } 772 }