1 /*
2 * R : A Computer Language for Statistical Data Analysis
3 * Copyright (C) 1995-2018 The R Core Team
4 * Copyright (C) 2003 The R Foundation
5 *
6 * This program is free software; you can redistribute it and/or modify
7 * it under the terms of the GNU General Public License as published by
8 * the Free Software Foundation; either version 2 of the License, or
9 * (at your option) any later version.
10 *
11 * This program is distributed in the hope that it will be useful,
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 * GNU General Public License for more details.
15 *
16 * You should have received a copy of the GNU General Public License
17 * along with this program; if not, a copy is available at
18 * https://www.R-project.org/Licenses/
19 */
20
21 #ifdef HAVE_CONFIG_H
22 #include <config.h>
23 #endif
24
25 #ifdef HAVE_LONG_DOUBLE
26 # define SQRTL sqrtl
27 #else
28 # define SQRTL sqrt
29 #endif
30
31 #include <Defn.h>
32 #include <Rmath.h>
33
34 #include "statsR.h"
35 #undef _
36 #ifdef ENABLE_NLS
37 #include <libintl.h>
38 #define _(String) dgettext ("stats", String)
39 #else
40 #define _(String) (String)
41 #endif
42
43 static SEXP corcov(SEXP x, SEXP y, SEXP na_method, SEXP kendall, Rboolean cor);
44
45
cor(SEXP x,SEXP y,SEXP na_method,SEXP kendall)46 SEXP cor(SEXP x, SEXP y, SEXP na_method, SEXP kendall)
47 {
48 return corcov(x, y, na_method, kendall, TRUE);
49 }
cov(SEXP x,SEXP y,SEXP na_method,SEXP kendall)50 SEXP cov(SEXP x, SEXP y, SEXP na_method, SEXP kendall)
51 {
52 return corcov(x, y, na_method, kendall, FALSE);
53 }
54
55
56
57 #define COV_SUM_UPDATE \
58 sum += xm * ym; \
59 if(cor) { \
60 xsd += xm * xm; \
61 ysd += ym * ym; \
62 }
63
64 #define ANS(I,J) ans[I + J * ncx]
65 #define CLAMP(X) (X >= 1. ? 1. : (X <= -1. ? -1. : X))
66
67 /* Note that "if (kendall)" and "if (cor)" are used inside a double for() loop;
68 which makes the code better readable -- and is hopefully dealt with
69 by a smartly optimizing compiler
70 */
71
72 /** Compute Cov(xx[], yy[]) or Cor(.,.) with n = length(xx)
73 */
74 #define COV_PAIRWISE_BODY \
75 LDOUBLE sum, xmean = 0., ymean = 0., xsd, ysd, xm, ym; \
76 int k, nobs, n1 = -1; /* -Wall initializing */ \
77 \
78 nobs = 0; \
79 if(!kendall) { \
80 xmean = ymean = 0.; \
81 for (k = 0 ; k < n ; k++) { \
82 if(!(ISNAN(xx[k]) || ISNAN(yy[k]))) { \
83 nobs ++; \
84 xmean += xx[k]; \
85 ymean += yy[k]; \
86 } \
87 } \
88 } else /*kendall*/ \
89 for (k = 0 ; k < n ; k++) \
90 if(!(ISNAN(xx[k]) || ISNAN(yy[k]))) \
91 nobs ++; \
92 \
93 if (nobs >= 2) { \
94 xsd = ysd = sum = 0.; \
95 if(!kendall) { \
96 xmean /= nobs; \
97 ymean /= nobs; \
98 n1 = nobs-1; \
99 } \
100 for(k=0; k < n; k++) { \
101 if(!(ISNAN(xx[k]) || ISNAN(yy[k]))) { \
102 if(!kendall) { \
103 xm = xx[k] - xmean; \
104 ym = yy[k] - ymean; \
105 \
106 COV_SUM_UPDATE \
107 } \
108 else { /* Kendall's tau */ \
109 for(n1=0 ; n1 < k ; n1++) \
110 if(!(ISNAN(xx[n1]) || ISNAN(yy[n1]))) { \
111 xm = sign(xx[k] - xx[n1]); \
112 ym = sign(yy[k] - yy[n1]); \
113 \
114 COV_SUM_UPDATE \
115 } \
116 } \
117 } \
118 } \
119 if (cor) { \
120 if(xsd == 0. || ysd == 0.) { \
121 *sd_0 = TRUE; \
122 sum = NA_REAL; \
123 } \
124 else { \
125 if(!kendall) { \
126 xsd /= n1; \
127 ysd /= n1; \
128 sum /= n1; \
129 } \
130 sum /= (SQRTL(xsd) * SQRTL(ysd)); \
131 sum = CLAMP(sum); \
132 } \
133 } \
134 else if(!kendall) \
135 sum /= n1; \
136 \
137 ANS(i,j) = (double) sum; \
138 } \
139 else \
140 ANS(i,j) = NA_REAL
141
142
cov_pairwise1(int n,int ncx,double * x,double * ans,Rboolean * sd_0,Rboolean cor,Rboolean kendall)143 static void cov_pairwise1(int n, int ncx, double *x,
144 double *ans, Rboolean *sd_0, Rboolean cor,
145 Rboolean kendall)
146 {
147 for (int i = 0 ; i < ncx ; i++) {
148 double *xx = &x[i * n];
149 for (int j = 0 ; j <= i ; j++) {
150 double *yy = &x[j * n];
151
152 COV_PAIRWISE_BODY;
153
154 ANS(j,i) = ANS(i,j);
155 }
156 }
157 }
158
cov_pairwise2(int n,int ncx,int ncy,double * x,double * y,double * ans,Rboolean * sd_0,Rboolean cor,Rboolean kendall)159 static void cov_pairwise2(int n, int ncx, int ncy, double *x, double *y,
160 double *ans, Rboolean *sd_0, Rboolean cor,
161 Rboolean kendall)
162 {
163 for (int i = 0 ; i < ncx ; i++) {
164 double *xx = &x[i * n];
165 for (int j = 0 ; j < ncy ; j++) {
166 double *yy = &y[j * n];
167
168 COV_PAIRWISE_BODY;
169 }
170 }
171 }
172 #undef COV_PAIRWISE_BODY
173
174
175 /* method = "complete" or "all.obs" (only difference: na_fail):
176 * -------- -------
177 */
178 #define COV_ini_0 \
179 LDOUBLE sum, tmp, xxm, yym; \
180 double *xx, *yy; \
181 R_xlen_t i, j, k, n1=-1/* -Wall */
182
183 #define COV_n_le_1(_n_,_k_) \
184 if (_n_ <= 1) {/* too many missing */ \
185 for (i = 0 ; i < ncx ; i++) \
186 for (j = 0 ; j < _k_ ; j++) \
187 ANS(i,j) = NA_REAL; \
188 return; \
189 }
190
191 #define COV_init(_ny_) \
192 COV_ini_0; int nobs; \
193 \
194 /* total number of complete observations */ \
195 nobs = 0; \
196 for(k = 0 ; k < n ; k++) { \
197 if (ind[k] != 0) nobs++; \
198 } \
199 COV_n_le_1(nobs, _ny_)
200
201 #define COV_ini_na(_ny_) \
202 COV_ini_0; \
203 COV_n_le_1(n, _ny_)
204
205
206 /* This uses two passes for better accuracy */
207 #define MEAN(_X_) \
208 /* variable means */ \
209 for (i = 0 ; i < nc##_X_ ; i++) { \
210 xx = &_X_[i * n]; \
211 sum = 0.; \
212 for (k = 0 ; k < n ; k++) \
213 if(ind[k] != 0) \
214 sum += xx[k]; \
215 tmp = sum / nobs; \
216 if(R_FINITE((double)tmp)) { \
217 sum = 0.; \
218 for (k = 0 ; k < n ; k++) \
219 if(ind[k] != 0) \
220 sum += (xx[k] - tmp); \
221 tmp = tmp + sum / nobs; \
222 } \
223 _X_##m [i] = (double)tmp; \
224 }
225
226 /* This uses two passes for better accuracy */
227 #define MEAN_(_X_,_HAS_NA_) \
228 /* variable means (has_na) */ \
229 for (i = 0 ; i < nc##_X_ ; i++) { \
230 if(_HAS_NA_[i]) \
231 tmp = NA_REAL; \
232 else { \
233 xx = &_X_[i * n]; \
234 sum = 0.; \
235 for (k = 0 ; k < n ; k++) \
236 sum += xx[k]; \
237 tmp = sum / n; \
238 if(R_FINITE((double)tmp)) { \
239 sum = 0.; \
240 for (k = 0 ; k < n ; k++) \
241 sum += (xx[k] - tmp); \
242 tmp = tmp + sum / n; \
243 } \
244 } \
245 _X_##m [i] = (double)tmp; \
246 }
247
248
249 static void
cov_complete1(int n,int ncx,double * x,double * xm,int * ind,double * ans,Rboolean * sd_0,Rboolean cor,Rboolean kendall)250 cov_complete1(int n, int ncx, double *x, double *xm,
251 int *ind, double *ans, Rboolean *sd_0, Rboolean cor,
252 Rboolean kendall)
253 {
254 COV_init(ncx);
255
256 if(!kendall) {
257 MEAN(x);/* -> xm[] */
258 n1 = nobs - 1;
259 }
260 for (i = 0 ; i < ncx ; i++) {
261 xx = &x[i * n];
262
263 if(!kendall) {
264 xxm = xm[i];
265 for (j = 0 ; j <= i ; j++) {
266 yy = &x[j * n];
267 yym = xm[j];
268 sum = 0.;
269 for (k = 0 ; k < n ; k++)
270 if (ind[k] != 0)
271 sum += (LDOUBLE)(xx[k] - xxm) * (yy[k] - yym);
272 ANS(j,i) = ANS(i,j) = (double)(sum / n1);
273 }
274 }
275 else { /* Kendall's tau */
276 for (j = 0 ; j <= i ; j++) {
277 yy = &x[j * n];
278 sum = 0.;
279 for (k = 0 ; k < n ; k++)
280 if (ind[k] != 0)
281 for (n1 = 0 ; n1 < n ; n1++)
282 if (ind[n1] != 0)
283 sum += sign(xx[k] - xx[n1])
284 * sign(yy[k] - yy[n1]);
285 ANS(j,i) = ANS(i,j) = (double)sum;
286 }
287 }
288 }
289
290 if (cor) {
291 for (i = 0 ; i < ncx ; i++)
292 xm[i] = sqrt(ANS(i,i));
293 for (i = 0 ; i < ncx ; i++) {
294 for (j = 0 ; j < i ; j++) {
295 if (xm[i] == 0 || xm[j] == 0) {
296 *sd_0 = TRUE;
297 ANS(j,i) = ANS(i,j) = NA_REAL;
298 }
299 else {
300 sum = ANS(i,j) / (xm[i] * xm[j]);
301 ANS(j,i) = ANS(i,j) = (double)CLAMP(sum);
302 }
303 }
304 ANS(i,i) = 1.0;
305 }
306 }
307 } /* cov_complete1 */
308
309 static void
cov_na_1(int n,int ncx,double * x,double * xm,int * has_na,double * ans,Rboolean * sd_0,Rboolean cor,Rboolean kendall)310 cov_na_1(int n, int ncx, double *x, double *xm,
311 int *has_na, double *ans, Rboolean *sd_0, Rboolean cor,
312 Rboolean kendall)
313 {
314
315 COV_ini_na(ncx);
316
317 if(!kendall) {
318 MEAN_(x, has_na);/* -> xm[] */
319 n1 = n - 1;
320 }
321 for (i = 0 ; i < ncx ; i++) {
322 if(has_na[i]) {
323 for (j = 0 ; j <= i ; j++)
324 ANS(j,i) = ANS(i,j) = NA_REAL;
325 }
326 else {
327 xx = &x[i * n];
328
329 if(!kendall) {
330 xxm = xm[i];
331 for (j = 0 ; j <= i ; j++)
332 if(has_na[j]) {
333 ANS(j,i) = ANS(i,j) = NA_REAL;
334 } else {
335 yy = &x[j * n];
336 yym = xm[j];
337 sum = 0.;
338 for (k = 0 ; k < n ; k++)
339 sum += (LDOUBLE)(xx[k] - xxm) * (yy[k] - yym);
340 ANS(j,i) = ANS(i,j) = (double)(sum / n1);
341 }
342 }
343 else { /* Kendall's tau */
344 for (j = 0 ; j <= i ; j++)
345 if(has_na[j]) {
346 ANS(j,i) = ANS(i,j) = NA_REAL;
347 } else {
348 yy = &x[j * n];
349 sum = 0.;
350 for (k = 0 ; k < n ; k++)
351 for (n1 = 0 ; n1 < n ; n1++)
352 sum += sign(xx[k] - xx[n1]) * sign(yy[k] - yy[n1]);
353 ANS(j,i) = ANS(i,j) = (double)sum;
354 }
355 }
356 }
357 }
358
359 if (cor) {
360 for (i = 0 ; i < ncx ; i++)
361 if(!has_na[i]) xm[i] = sqrt(ANS(i,i));
362 for (i = 0 ; i < ncx ; i++) {
363 if(!has_na[i]) for (j = 0 ; j < i ; j++) {
364 if (xm[i] == 0 || xm[j] == 0) {
365 *sd_0 = TRUE;
366 ANS(j,i) = ANS(i,j) = NA_REAL;
367 }
368 else {
369 sum = ANS(i,j) / (xm[i] * xm[j]);
370 ANS(j,i) = ANS(i,j) = (double)CLAMP(sum);
371 }
372 }
373 ANS(i,i) = 1.0;
374 }
375 }
376 } /* cov_na_1() */
377
378 static void
cov_complete2(int n,int ncx,int ncy,double * x,double * y,double * xm,double * ym,int * ind,double * ans,Rboolean * sd_0,Rboolean cor,Rboolean kendall)379 cov_complete2(int n, int ncx, int ncy, double *x, double *y,
380 double *xm, double *ym, int *ind,
381 double *ans, Rboolean *sd_0, Rboolean cor, Rboolean kendall)
382 {
383 COV_init(ncy);
384
385 if(!kendall) {
386 MEAN(x);/* -> xm[] */
387 MEAN(y);/* -> ym[] */
388 n1 = nobs - 1;
389 }
390 for (i = 0 ; i < ncx ; i++) {
391 xx = &x[i * n];
392 if(!kendall) {
393 xxm = xm[i];
394 for (j = 0 ; j < ncy ; j++) {
395 yy = &y[j * n];
396 yym = ym[j];
397 sum = 0.;
398 for (k = 0 ; k < n ; k++)
399 if (ind[k] != 0)
400 sum += (LDOUBLE)(xx[k] - xxm) * (yy[k] - yym);
401 ANS(i,j) = (double)(sum / n1);
402 }
403 }
404 else { /* Kendall's tau */
405 for (j = 0 ; j < ncy ; j++) {
406 yy = &y[j * n];
407 sum = 0.;
408 for (k = 0 ; k < n ; k++)
409 if (ind[k] != 0)
410 for (n1 = 0 ; n1 < n ; n1++)
411 if (ind[n1] != 0)
412 sum += sign(xx[k] - xx[n1])
413 * sign(yy[k] - yy[n1]);
414 ANS(i,j) = (double)sum;
415 }
416 }
417 }
418
419 if (cor) {
420
421 #define COV_SDEV(_X_) \
422 for (i = 0 ; i < nc##_X_ ; i++) { /* Var(X[i]) */ \
423 xx = &_X_[i * n]; \
424 sum = 0.; \
425 if(!kendall) { \
426 xxm = _X_##m [i]; \
427 for (k = 0 ; k < n ; k++) \
428 if (ind[k] != 0) \
429 sum += (LDOUBLE)(xx[k] - xxm) * (xx[k] - xxm); \
430 sum /= n1; \
431 } \
432 else { /* Kendall's tau */ \
433 for (k = 0 ; k < n ; k++) \
434 if (ind[k] != 0) \
435 for (n1 = 0 ; n1 < n ; n1++) \
436 if (ind[n1] != 0 && xx[k] != xx[n1]) \
437 sum ++; /* = sign(. - .)^2 */ \
438 } \
439 _X_##m [i] = (double)SQRTL(sum); \
440 }
441
442 COV_SDEV(x); /* -> xm[.] */
443 COV_SDEV(y); /* -> ym[.] */
444
445 for (i = 0 ; i < ncx ; i++)
446 for (j = 0 ; j < ncy ; j++)
447 if (xm[i] == 0. || ym[j] == 0.) {
448 *sd_0 = TRUE;
449 ANS(i,j) = NA_REAL;
450 }
451 else {
452 ANS(i,j) /= (xm[i] * ym[j]);
453 ANS(i,j) = CLAMP(ANS(i,j));
454 }
455 }/* cor */
456
457 }/* cov_complete2 */
458 #undef COV_SDEV
459
460 static void
cov_na_2(int n,int ncx,int ncy,double * x,double * y,double * xm,double * ym,int * has_na_x,int * has_na_y,double * ans,Rboolean * sd_0,Rboolean cor,Rboolean kendall)461 cov_na_2(int n, int ncx, int ncy, double *x, double *y,
462 double *xm, double *ym, int *has_na_x, int *has_na_y,
463 double *ans, Rboolean *sd_0, Rboolean cor, Rboolean kendall)
464 {
465 COV_ini_na(ncy);
466
467 if(!kendall) {
468 MEAN_(x, has_na_x);/* -> xm[] */
469 MEAN_(y, has_na_y);/* -> ym[] */
470 n1 = n - 1;
471 }
472 for (i = 0 ; i < ncx ; i++) {
473 if(has_na_x[i]) {
474 for (j = 0 ; j < ncy; j++)
475 ANS(i,j) = NA_REAL;
476 }
477 else {
478 xx = &x[i * n];
479 if(!kendall) {
480 xxm = xm[i];
481 for (j = 0 ; j < ncy ; j++)
482 if(has_na_y[j]) {
483 ANS(i,j) = NA_REAL;
484 } else {
485 yy = &y[j * n];
486 yym = ym[j];
487 sum = 0.;
488 for (k = 0 ; k < n ; k++)
489 sum += (LDOUBLE)(xx[k] - xxm) * (yy[k] - yym);
490 ANS(i,j) = (double)(sum / n1);
491 }
492 }
493 else { /* Kendall's tau */
494 for (j = 0 ; j < ncy ; j++)
495 if(has_na_y[j]) {
496 ANS(i,j) = NA_REAL;
497 } else {
498 yy = &y[j * n];
499 sum = 0.;
500 for (k = 0 ; k < n ; k++)
501 for (n1 = 0 ; n1 < n ; n1++)
502 sum += sign(xx[k] - xx[n1]) * sign(yy[k] - yy[n1]);
503 ANS(i,j) = (double)sum;
504 }
505 }
506 }
507 }
508
509 if (cor) {
510
511 #define COV_SDEV(_X_) \
512 for (i = 0 ; i < nc##_X_ ; i++) \
513 if(!has_na_##_X_[i]) { /* Var(X[j]) */ \
514 xx = &_X_[i * n]; \
515 sum = 0.; \
516 if(!kendall) { \
517 xxm = _X_##m [i]; \
518 for (k = 0 ; k < n ; k++) \
519 sum += (LDOUBLE)(xx[k] - xxm) * (xx[k] - xxm); \
520 sum /= n1; \
521 } \
522 else { /* Kendall's tau */ \
523 for (k = 0 ; k < n ; k++) \
524 for (n1 = 0 ; n1 < n ; n1++) \
525 if (xx[k] != xx[n1]) \
526 sum ++; /* = sign(. - .)^2 */ \
527 } \
528 _X_##m [i] = (double) SQRTL(sum); \
529 }
530
531 COV_SDEV(x); /* -> xm[.] */
532 COV_SDEV(y); /* -> ym[.] */
533
534 for (i = 0 ; i < ncx ; i++)
535 if(!has_na_x[i]) {
536 for (j = 0 ; j < ncy ; j++)
537 if(!has_na_y[j]) {
538 if (xm[i] == 0. || ym[j] == 0.) {
539 *sd_0 = TRUE;
540 ANS(i,j) = NA_REAL;
541 }
542 else {
543 ANS(i,j) /= (xm[i] * ym[j]);
544 ANS(i,j) = CLAMP(ANS(i,j));
545 }
546 }
547 }
548 }/* cor */
549
550 }/* cov_na_2 */
551
552 #undef ANS
553 #undef COV_init
554 #undef MEAN
555 #undef MEAN_
556 #undef COV_SDEV
557
558 /* complete[12]() returns indicator vector ind[] of complete.cases(), or
559 * -------------- if(na_fail) signals error if any NA/NaN is encountered
560 */
561
562 /* This might look slightly inefficient, but it is designed to
563 * optimise paging in virtual memory systems ...
564 * (or at least that's my story, and I'm sticking to it.)
565 */
566 #define NA_LOOP \
567 for (i = 0 ; i < n ; i++) \
568 if (ISNAN(z[i])) { \
569 if (na_fail) error(_("missing observations in cov/cor"));\
570 else ind[i] = 0; \
571 }
572
573 #define COMPLETE_1 \
574 double *z; \
575 int i, j; \
576 for (i = 0 ; i < n ; i++) \
577 ind[i] = 1; \
578 for (j = 0 ; j < ncx ; j++) { \
579 z = &x[j * n]; \
580 NA_LOOP \
581 }
582
complete1(int n,int ncx,double * x,int * ind,Rboolean na_fail)583 static void complete1(int n, int ncx, double *x, int *ind, Rboolean na_fail)
584 {
585 COMPLETE_1
586 }
587
588 static void
complete2(int n,int ncx,int ncy,double * x,double * y,int * ind,Rboolean na_fail)589 complete2(int n, int ncx, int ncy, double *x, double *y, int *ind, Rboolean na_fail)
590 {
591 COMPLETE_1
592
593 for(j = 0 ; j < ncy ; j++) {
594 z = &y[j * n];
595 NA_LOOP
596 }
597 }
598
599 #define HAS_NA_1(_X_,_HAS_NA_) \
600 for (j = 0 ; j < nc##_X_ ; j++) { \
601 z = &_X_[j * n]; \
602 _HAS_NA_[j] = 0; \
603 for (i = 0 ; i < n ; i++) \
604 if (ISNAN(z[i])) { \
605 _HAS_NA_[j] = 1; break; \
606 } \
607 }
608
609
find_na_1(int n,int ncx,double * x,int * has_na)610 static void find_na_1(int n, int ncx, double *x, int *has_na)
611 {
612 double *z;
613 int i, j;
614 HAS_NA_1(x, has_na)
615 }
616
617 static void
find_na_2(int n,int ncx,int ncy,double * x,double * y,int * has_na_x,int * has_na_y)618 find_na_2(int n, int ncx, int ncy, double *x, double *y, int *has_na_x, int *has_na_y)
619 {
620 double *z;
621 int i, j;
622 HAS_NA_1(x, has_na_x)
623 HAS_NA_1(y, has_na_y)
624 }
625
626 #undef NA_LOOP
627 #undef COMPLETE_1
628 #undef HAS_NA_1
629
630 /* co[vr](x, y, use =
631 { 1, 2, 3, 4, 5 }
632 "all.obs", "complete.obs", "pairwise.complete", "everything", "na.or.complete"
633 kendall = TRUE/FALSE)
634 */
corcov(SEXP x,SEXP y,SEXP na_method,SEXP skendall,Rboolean cor)635 static SEXP corcov(SEXP x, SEXP y, SEXP na_method, SEXP skendall, Rboolean cor)
636 {
637 SEXP ans, xm, ym, ind;
638 Rboolean ansmat, kendall, pair, na_fail, everything, sd_0, empty_err;
639 int i, method, n, ncx, ncy, nprotect = 2;
640
641 #define DEFUNCT_VAR_FACTOR
642 #ifdef DEFUNCT_VAR_FACTOR
643 # define VAR_FACTOR_MSG "Calling var(x) on a factor x is defunct.\n Use something like 'all(duplicated(x)[-1L])' to test for a constant vector."
644 #else
645 # define VAR_FACTOR_MSG "Calling var(x) on a factor x is deprecated and will become an error.\n Use something like 'all(duplicated(x)[-1L])' to test for a constant vector."
646 #endif
647
648 /* Arg.1: x */
649 if(isNull(x)) /* never allowed */
650 error(_("'x' is NULL"));
651 if(isFactor(x))
652 #ifdef DEFUNCT_VAR_FACTOR
653 error(_(VAR_FACTOR_MSG));
654 #else
655 warning(_(VAR_FACTOR_MSG));
656 #endif
657
658 /* length check of x -- only if(empty_err) --> below */
659 x = PROTECT(coerceVector(x, REALSXP));
660 if ((ansmat = isMatrix(x))) {
661 n = nrows(x);
662 ncx = ncols(x);
663 }
664 else {
665 n = length(x);
666 ncx = 1;
667 }
668 /* Arg.2: y */
669 if (isNull(y)) {/* y = x : var() */
670 ncy = ncx;
671 } else {
672 if(isFactor(y))
673 #ifdef DEFUNCT_VAR_FACTOR
674 error(_(VAR_FACTOR_MSG));
675 #else
676 warning(_(VAR_FACTOR_MSG));
677 #endif
678 y = PROTECT(coerceVector(y, REALSXP));
679 nprotect++;
680 if (isMatrix(y)) {
681 if (nrows(y) != n)
682 error(_("incompatible dimensions"));
683 ncy = ncols(y);
684 ansmat = TRUE;
685 }
686 else {
687 if (length(y) != n)
688 error(_("incompatible dimensions"));
689 ncy = 1;
690 }
691 }
692 /* Arg.3: method */
693 method = asInteger(na_method);
694
695 /* Arg.4: kendall */
696 kendall = asLogical(skendall);
697
698 /* "default: complete" (easier for -Wall) */
699 na_fail = FALSE; everything = FALSE; empty_err = TRUE;
700 pair = FALSE;
701 switch(method) {
702 case 1: /* use all : no NAs */
703 na_fail = TRUE;
704 break;
705 case 2: /* complete */
706 /* did na.omit in R */
707 if (!LENGTH(x)) error(_("no complete element pairs"));
708 break;
709 case 3: /* pairwise.complete */
710 pair = TRUE;
711 break;
712 case 4: /* "everything": NAs are propagated */
713 everything = TRUE;
714 empty_err = FALSE;
715 break;
716 case 5: /* "na.or.complete": NAs are propagated */
717 empty_err = FALSE;
718 break;
719 default:
720 error(_("invalid 'use' (computational method)"));
721 }
722 if (empty_err && !LENGTH(x))
723 error(_("'x' is empty"));
724
725 if (ansmat) PROTECT(ans = allocMatrix(REALSXP, ncx, ncy));
726 else PROTECT(ans = allocVector(REALSXP, ncx * ncy));
727 sd_0 = FALSE;
728 if (isNull(y)) {
729 if (everything) { /* NA's are propagated */
730 PROTECT(xm = allocVector(REALSXP, ncx));
731 PROTECT(ind = allocVector(LGLSXP, ncx));
732 find_na_1(n, ncx, REAL(x), /* --> has_na[] = */ LOGICAL(ind));
733 cov_na_1 (n, ncx, REAL(x), REAL(xm), LOGICAL(ind), REAL(ans), &sd_0, cor, kendall);
734
735 UNPROTECT(2);
736 }
737 else if (!pair) { /* all | complete "var" */
738 PROTECT(xm = allocVector(REALSXP, ncx));
739 PROTECT(ind = allocVector(INTSXP, n));
740 complete1(n, ncx, REAL(x), INTEGER(ind), na_fail);
741 cov_complete1(n, ncx, REAL(x), REAL(xm),
742 INTEGER(ind), REAL(ans), &sd_0, cor, kendall);
743 if(empty_err) {
744 Rboolean indany = FALSE;
745 for(i = 0; i < n; i++) {
746 if(INTEGER(ind)[i] == 1) { indany = TRUE; break; }
747 }
748 if(!indany) error(_("no complete element pairs"));
749 }
750 UNPROTECT(2);
751 }
752 else { /* pairwise "var" */
753 cov_pairwise1(n, ncx, REAL(x), REAL(ans), &sd_0, cor, kendall);
754 }
755 }
756 else { /* Co[vr] (x, y) */
757 if (everything) {
758 SEXP has_na_y;
759 PROTECT(xm = allocVector(REALSXP, ncx));
760 PROTECT(ym = allocVector(REALSXP, ncy));
761 PROTECT(ind = allocVector(LGLSXP, ncx));
762 PROTECT(has_na_y = allocVector(LGLSXP, ncy));
763
764 find_na_2(n, ncx, ncy, REAL(x), REAL(y), INTEGER(ind), INTEGER(has_na_y));
765 cov_na_2 (n, ncx, ncy, REAL(x), REAL(y), REAL(xm), REAL(ym),
766 INTEGER(ind), INTEGER(has_na_y), REAL(ans), &sd_0, cor, kendall);
767 UNPROTECT(4);
768 }
769 else if (!pair) { /* all | complete */
770 PROTECT(xm = allocVector(REALSXP, ncx));
771 PROTECT(ym = allocVector(REALSXP, ncy));
772 PROTECT(ind = allocVector(INTSXP, n));
773 complete2(n, ncx, ncy, REAL(x), REAL(y), INTEGER(ind), na_fail);
774 cov_complete2(n, ncx, ncy, REAL(x), REAL(y), REAL(xm), REAL(ym),
775 INTEGER(ind), REAL(ans), &sd_0, cor, kendall);
776 if(empty_err) {
777 Rboolean indany = FALSE;
778 for(i = 0; i < n; i++) {
779 if(INTEGER(ind)[i] == 1) { indany = TRUE; break; }
780 }
781 if(!indany) error(_("no complete element pairs"));
782 }
783 UNPROTECT(3);
784 }
785 else { /* pairwise */
786 cov_pairwise2(n, ncx, ncy, REAL(x), REAL(y), REAL(ans),
787 &sd_0, cor, kendall);
788 }
789 }
790 if (ansmat) { /* set dimnames() when applicable */
791 if (isNull(y)) {
792 x = getAttrib(x, R_DimNamesSymbol);
793 if (!isNull(x) && !isNull(VECTOR_ELT(x, 1))) {
794 PROTECT(ind = allocVector(VECSXP, 2));
795 SET_VECTOR_ELT(ind, 0, duplicate(VECTOR_ELT(x, 1)));
796 SET_VECTOR_ELT(ind, 1, duplicate(VECTOR_ELT(x, 1)));
797 setAttrib(ans, R_DimNamesSymbol, ind);
798 UNPROTECT(1);
799 }
800 }
801 else {
802 x = getAttrib(x, R_DimNamesSymbol);
803 y = getAttrib(y, R_DimNamesSymbol);
804 if ((length(x) >= 2 && !isNull(VECTOR_ELT(x, 1))) ||
805 (length(y) >= 2 && !isNull(VECTOR_ELT(y, 1)))) {
806 PROTECT(ind = allocVector(VECSXP, 2));
807 if (length(x) >= 2 && !isNull(VECTOR_ELT(x, 1)))
808 SET_VECTOR_ELT(ind, 0, duplicate(VECTOR_ELT(x, 1)));
809 if (length(y) >= 2 && !isNull(VECTOR_ELT(y, 1)))
810 SET_VECTOR_ELT(ind, 1, duplicate(VECTOR_ELT(y, 1)));
811 setAttrib(ans, R_DimNamesSymbol, ind);
812 UNPROTECT(1);
813 }
814 }
815 }
816 if(sd_0)/* only in cor() */
817 warning(_("the standard deviation is zero"));
818 UNPROTECT(nprotect);
819 return ans;
820 }
821