R-  3.3.1
About: R (also known as "GNU S") is a system for statistical computation and graphics.
  Fossies Dox: R-3.3.1.tar.gz  ("inofficial" and yet experimental doxygen-generated source code documentation)  

cov.c
Go to the documentation of this file.
1 /*
2  * R : A Computer Language for Statistical Data Analysis
3  * Copyright (C) 1995-2015 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 
46 SEXP cor(SEXP x, SEXP y, SEXP na_method, SEXP kendall)
47 {
48  return corcov(x, y, na_method, kendall, TRUE);
49 }
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 
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 
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 
153 
154  ANS(j,i) = ANS(i,j);
155  }
156  }
157 }
158 
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 
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  int 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
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 += (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
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 += (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
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 += (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 += (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
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 += (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 += (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 
583 static void complete1(int n, int ncx, double *x, int *ind, Rboolean na_fail)
584 {
585  COMPLETE_1
586 }
587 
588 static void
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 
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
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 */
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  /* Arg.1: x */
642  if(isNull(x)) /* never allowed */
643  error(_("'x' is NULL"));
644 #ifdef _R_in_2017_
645  if(isFactor(x)) error(_("'x' is a factor"));
646 #else
647 # 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."
648  if(isFactor(x)) warning(_(VAR_FACTOR_MSG));
649 #endif
650  /* length check of x -- only if(empty_err) --> below */
651  x = PROTECT(coerceVector(x, REALSXP));
652  if ((ansmat = isMatrix(x))) {
653  n = nrows(x);
654  ncx = ncols(x);
655  }
656  else {
657  n = length(x);
658  ncx = 1;
659  }
660  /* Arg.2: y */
661  if (isNull(y)) {/* y = x : var() */
662  ncy = ncx;
663  } else {
664 #ifdef _R_in_2017_
665  if(isFactor(y)) error(_("'y' is a factor"));
666 #else
668 #endif
670  nprotect++;
671  if (isMatrix(y)) {
672  if (nrows(y) != n)
673  error(_("incompatible dimensions"));
674  ncy = ncols(y);
675  ansmat = TRUE;
676  }
677  else {
678  if (length(y) != n)
679  error(_("incompatible dimensions"));
680  ncy = 1;
681  }
682  }
683  /* Arg.3: method */
684  method = asInteger(na_method);
685 
686  /* Arg.4: kendall */
687  kendall = asLogical(skendall);
688 
689  /* "default: complete" (easier for -Wall) */
690  na_fail = FALSE; everything = FALSE; empty_err = TRUE;
691  pair = FALSE;
692  switch(method) {
693  case 1: /* use all : no NAs */
694  na_fail = TRUE;
695  break;
696  case 2: /* complete */
697  /* did na.omit in R */
698  if (!LENGTH(x)) error(_("no complete element pairs"));
699  break;
700  case 3: /* pairwise.complete */
701  pair = TRUE;
702  break;
703  case 4: /* "everything": NAs are propagated */
704  everything = TRUE;
705  empty_err = FALSE;
706  break;
707  case 5: /* "na.or.complete": NAs are propagated */
708  empty_err = FALSE;
709  break;
710  default:
711  error(_("invalid 'use' (computational method)"));
712  }
713  if (empty_err && !LENGTH(x))
714  error(_("'x' is empty"));
715 
716  if (ansmat) PROTECT(ans = allocMatrix(REALSXP, ncx, ncy));
717  else PROTECT(ans = allocVector(REALSXP, ncx * ncy));
718  sd_0 = FALSE;
719  if (isNull(y)) {
720  if (everything) { /* NA's are propagated */
721  PROTECT(xm = allocVector(REALSXP, ncx));
722  PROTECT(ind = allocVector(LGLSXP, ncx));
723  find_na_1(n, ncx, REAL(x), /* --> has_na[] = */ LOGICAL(ind));
724  cov_na_1 (n, ncx, REAL(x), REAL(xm), LOGICAL(ind), REAL(ans), &sd_0, cor, kendall);
725 
726  UNPROTECT(2);
727  }
728  else if (!pair) { /* all | complete "var" */
729  PROTECT(xm = allocVector(REALSXP, ncx));
730  PROTECT(ind = allocVector(INTSXP, n));
731  complete1(n, ncx, REAL(x), INTEGER(ind), na_fail);
732  cov_complete1(n, ncx, REAL(x), REAL(xm),
733  INTEGER(ind), REAL(ans), &sd_0, cor, kendall);
734  if(empty_err) {
735  Rboolean indany = FALSE;
736  for(i = 0; i < n; i++) {
737  if(INTEGER(ind)[i] == 1) { indany = TRUE; break; }
738  }
739  if(!indany) error(_("no complete element pairs"));
740  }
741  UNPROTECT(2);
742  }
743  else { /* pairwise "var" */
744  cov_pairwise1(n, ncx, REAL(x), REAL(ans), &sd_0, cor, kendall);
745  }
746  }
747  else { /* Co[vr] (x, y) */
748  if (everything) {
749  SEXP has_na_y;
750  PROTECT(xm = allocVector(REALSXP, ncx));
751  PROTECT(ym = allocVector(REALSXP, ncy));
752  PROTECT(ind = allocVector(LGLSXP, ncx));
753  PROTECT(has_na_y = allocVector(LGLSXP, ncy));
754 
755  find_na_2(n, ncx, ncy, REAL(x), REAL(y), INTEGER(ind), INTEGER(has_na_y));
756  cov_na_2 (n, ncx, ncy, REAL(x), REAL(y), REAL(xm), REAL(ym),
757  INTEGER(ind), INTEGER(has_na_y), REAL(ans), &sd_0, cor, kendall);
758  UNPROTECT(4);
759  }
760  else if (!pair) { /* all | complete */
761  PROTECT(xm = allocVector(REALSXP, ncx));
762  PROTECT(ym = allocVector(REALSXP, ncy));
763  PROTECT(ind = allocVector(INTSXP, n));
764  complete2(n, ncx, ncy, REAL(x), REAL(y), INTEGER(ind), na_fail);
765  cov_complete2(n, ncx, ncy, REAL(x), REAL(y), REAL(xm), REAL(ym),
766  INTEGER(ind), REAL(ans), &sd_0, cor, kendall);
767  if(empty_err) {
768  Rboolean indany = FALSE;
769  for(i = 0; i < n; i++) {
770  if(INTEGER(ind)[i] == 1) { indany = TRUE; break; }
771  }
772  if(!indany) error(_("no complete element pairs"));
773  }
774  UNPROTECT(3);
775  }
776  else { /* pairwise */
777  cov_pairwise2(n, ncx, ncy, REAL(x), REAL(y), REAL(ans),
778  &sd_0, cor, kendall);
779  }
780  }
781  if (ansmat) { /* set dimnames() when applicable */
782  if (isNull(y)) {
783  x = getAttrib(x, R_DimNamesSymbol);
784  if (!isNull(x) && !isNull(VECTOR_ELT(x, 1))) {
785  PROTECT(ind = allocVector(VECSXP, 2));
786  SET_VECTOR_ELT(ind, 0, duplicate(VECTOR_ELT(x, 1)));
787  SET_VECTOR_ELT(ind, 1, duplicate(VECTOR_ELT(x, 1)));
788  setAttrib(ans, R_DimNamesSymbol, ind);
789  UNPROTECT(1);
790  }
791  }
792  else {
793  x = getAttrib(x, R_DimNamesSymbol);
795  if ((length(x) >= 2 && !isNull(VECTOR_ELT(x, 1))) ||
796  (length(y) >= 2 && !isNull(VECTOR_ELT(y, 1)))) {
797  PROTECT(ind = allocVector(VECSXP, 2));
798  if (length(x) >= 2 && !isNull(VECTOR_ELT(x, 1)))
799  SET_VECTOR_ELT(ind, 0, duplicate(VECTOR_ELT(x, 1)));
800  if (length(y) >= 2 && !isNull(VECTOR_ELT(y, 1)))
801  SET_VECTOR_ELT(ind, 1, duplicate(VECTOR_ELT(y, 1)));
802  setAttrib(ans, R_DimNamesSymbol, ind);
803  UNPROTECT(1);
804  }
805  }
806  }
807  if(sd_0)/* only in cor() */
808  warning(_("the standard deviation is zero"));
809  UNPROTECT(nprotect);
810  return ans;
811 }
#define VECSXP
Definition: Rinternals.h:131
#define NA_REAL
Definition: Arith.h:54
#define UNPROTECT(n)
Definition: Rinternals.h:662
LibExtern SEXP R_DimNamesSymbol
Definition: Rinternals.h:703
double *() REAL(SEXP x)
Definition: memory.c:3442
Rboolean
Definition: Boolean.h:31
#define REALSXP
Definition: Rinternals.h:124
SEXP SET_VECTOR_ELT(SEXP x, R_xlen_t i, SEXP v)
Definition: memory.c:3478
#define MEAN_(_X_, _HAS_NA_)
Definition: cov.c:227
SEXP coerceVector(SEXP v, SEXPTYPE type)
Definition: coerce.c:1148
#define HAS_NA_1(_X_, _HAS_NA_)
Definition: cov.c:599
void warning(const char *format,...)
Definition: errors.c:254
SEXP() VECTOR_ELT(SEXP x, R_xlen_t i)
Definition: memory.c:3410
j
Definition: qsort-body.c:133
#define COV_init(_ny_)
Definition: cov.c:191
#define y(k)
#define _(String)
Definition: cov.c:40
INLINE_FUN R_len_t length(SEXP s)
Definition: Rinlinedfuns.h:124
INLINE_FUN Rboolean isMatrix(SEXP s)
Definition: Rinlinedfuns.h:502
static void n
Definition: printvector.c:263
void error(const char *format,...)
Definition: errors.c:750
int *() LOGICAL(SEXP x)
Definition: memory.c:3420
SEXP setAttrib(SEXP vec, SEXP name, SEXP val)
Definition: attrib.c:214
SEXP cor(SEXP x, SEXP y, SEXP na_method, SEXP kendall)
Definition: cov.c:46
int asLogical(SEXP x)
Definition: coerce.c:1614
int nrows(SEXP s)
Definition: util.c:76
int ncols(SEXP s)
Definition: util.c:92
#define COMPLETE_1
Definition: cov.c:573
#define MEAN(_X_)
Definition: cov.c:207
#define COV_SDEV(_X_)
#define isNull
Definition: Rinternals.h:1192
#define ANS(I, J)
Definition: cov.c:64
int asInteger(SEXP x)
Definition: coerce.c:1643
#define COV_ini_na(_ny_)
Definition: cov.c:201
INTt k
Definition: qsort-body.c:48
Definition: Boolean.h:31
int() LENGTH(SEXP x)
Definition: memory.c:3388
#define NA_LOOP
Definition: cov.c:566
#define COV_PAIRWISE_BODY
Definition: cov.c:74
Definition: Boolean.h:31
struct SEXPREC * SEXP
Definition: Rinternals.h:446
#define PROTECT(s)
Definition: Rinternals.h:661
#define LGLSXP
Definition: Rinternals.h:121
SEXP getAttrib(SEXP vec, SEXP name)
Definition: attrib.c:151
double sign(double x)
Definition: sign.c:35
SEXP cov(SEXP x, SEXP y, SEXP na_method, SEXP kendall)
Definition: cov.c:50
INLINE_FUN Rboolean isFactor(SEXP s)
Definition: Rinlinedfuns.h:539
SEXP duplicate(SEXP s)
Definition: duplicate.c:128
#define CLAMP(X)
Definition: cov.c:65
INLINE_FUN SEXP allocVector(SEXPTYPE type, R_xlen_t length)
Definition: Rinlinedfuns.h:194
#define INTSXP
Definition: Rinternals.h:123
SEXP allocMatrix(SEXPTYPE mode, int nrow, int ncol)
Definition: array.c:202
#define VAR_FACTOR_MSG
int *() INTEGER(SEXP x)
Definition: memory.c:3428
int i
Definition: uuidgen.c:9