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