R-  3.2.3
About: R (also known as "GNU S") is a system for statistical computation and graphics.
  Fossies Dox: R-3.2.3.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 
66 /* Note that "if (kendall)" and "if (cor)" are used inside a double for() loop;
67  which makes the code better readable -- and is hopefully dealt with
68  by a smartly optimizing compiler
69 */
70 
73 #define COV_PAIRWISE_BODY \
74  LDOUBLE sum, xmean = 0., ymean = 0., xsd, ysd, xm, ym; \
75  int k, nobs, n1 = -1; /* -Wall initializing */ \
76  \
77  nobs = 0; \
78  if(!kendall) { \
79  xmean = ymean = 0.; \
80  for (k = 0 ; k < n ; k++) { \
81  if(!(ISNAN(xx[k]) || ISNAN(yy[k]))) { \
82  nobs ++; \
83  xmean += xx[k]; \
84  ymean += yy[k]; \
85  } \
86  } \
87  } else /*kendall*/ \
88  for (k = 0 ; k < n ; k++) \
89  if(!(ISNAN(xx[k]) || ISNAN(yy[k]))) \
90  nobs ++; \
91  \
92  if (nobs >= 2) { \
93  xsd = ysd = sum = 0.; \
94  if(!kendall) { \
95  xmean /= nobs; \
96  ymean /= nobs; \
97  n1 = nobs-1; \
98  } \
99  for(k=0; k < n; k++) { \
100  if(!(ISNAN(xx[k]) || ISNAN(yy[k]))) { \
101  if(!kendall) { \
102  xm = xx[k] - xmean; \
103  ym = yy[k] - ymean; \
104  \
105  COV_SUM_UPDATE \
106  } \
107  else { /* Kendall's tau */ \
108  for(n1=0 ; n1 < k ; n1++) \
109  if(!(ISNAN(xx[n1]) || ISNAN(yy[n1]))) { \
110  xm = sign(xx[k] - xx[n1]); \
111  ym = sign(yy[k] - yy[n1]); \
112  \
113  COV_SUM_UPDATE \
114  } \
115  } \
116  } \
117  } \
118  if (cor) { \
119  if(xsd == 0. || ysd == 0.) { \
120  *sd_0 = TRUE; \
121  sum = NA_REAL; \
122  } \
123  else { \
124  if(!kendall) { \
125  xsd /= n1; \
126  ysd /= n1; \
127  sum /= n1; \
128  } \
129  sum /= (SQRTL(xsd) * SQRTL(ysd)); \
130  if(sum > 1.) sum = 1.; \
131  } \
132  } \
133  else if(!kendall) \
134  sum /= n1; \
135  \
136  ANS(i,j) = (double) sum; \
137  } \
138  else \
139  ANS(i,j) = NA_REAL
140 
141 
142 static void cov_pairwise1(int n, int ncx, double *x,
143  double *ans, Rboolean *sd_0, Rboolean cor,
144  Rboolean kendall)
145 {
146  for (int i = 0 ; i < ncx ; i++) {
147  double *xx = &x[i * n];
148  for (int j = 0 ; j <= i ; j++) {
149  double *yy = &x[j * n];
150 
152 
153  ANS(j,i) = ANS(i,j);
154  }
155  }
156 }
157 
158 static void cov_pairwise2(int n, int ncx, int ncy, double *x, double *y,
159  double *ans, Rboolean *sd_0, Rboolean cor,
160  Rboolean kendall)
161 {
162  for (int i = 0 ; i < ncx ; i++) {
163  double *xx = &x[i * n];
164  for (int j = 0 ; j < ncy ; j++) {
165  double *yy = &y[j * n];
166 
168  }
169  }
170 }
171 #undef COV_PAIRWISE_BODY
172 
173 
174 /* method = "complete" or "all.obs" (only difference: na_fail):
175  * -------- -------
176 */
177 #define COV_ini_0 \
178  LDOUBLE sum, tmp, xxm, yym; \
179  double *xx, *yy; \
180  int i, j, k, n1=-1/* -Wall */
181 
182 #define COV_n_le_1(_n_,_k_) \
183  if (_n_ <= 1) {/* too many missing */ \
184  for (i = 0 ; i < ncx ; i++) \
185  for (j = 0 ; j < _k_ ; j++) \
186  ANS(i,j) = NA_REAL; \
187  return; \
188  }
189 
190 #define COV_init(_ny_) \
191  COV_ini_0; int nobs; \
192  \
193  /* total number of complete observations */ \
194  nobs = 0; \
195  for(k = 0 ; k < n ; k++) { \
196  if (ind[k] != 0) nobs++; \
197  } \
198  COV_n_le_1(nobs, _ny_)
199 
200 #define COV_ini_na(_ny_) \
201  COV_ini_0; \
202  COV_n_le_1(n, _ny_)
203 
204 
205 /* This uses two passes for better accuracy */
206 #define MEAN(_X_) \
207  /* variable means */ \
208  for (i = 0 ; i < nc##_X_ ; i++) { \
209  xx = &_X_[i * n]; \
210  sum = 0.; \
211  for (k = 0 ; k < n ; k++) \
212  if(ind[k] != 0) \
213  sum += xx[k]; \
214  tmp = sum / nobs; \
215  if(R_FINITE((double)tmp)) { \
216  sum = 0.; \
217  for (k = 0 ; k < n ; k++) \
218  if(ind[k] != 0) \
219  sum += (xx[k] - tmp); \
220  tmp = tmp + sum / nobs; \
221  } \
222  _X_##m [i] = (double)tmp; \
223  }
224 
225 /* This uses two passes for better accuracy */
226 #define MEAN_(_X_,_HAS_NA_) \
227  /* variable means (has_na) */ \
228  for (i = 0 ; i < nc##_X_ ; i++) { \
229  if(_HAS_NA_[i]) \
230  tmp = NA_REAL; \
231  else { \
232  xx = &_X_[i * n]; \
233  sum = 0.; \
234  for (k = 0 ; k < n ; k++) \
235  sum += xx[k]; \
236  tmp = sum / n; \
237  if(R_FINITE((double)tmp)) { \
238  sum = 0.; \
239  for (k = 0 ; k < n ; k++) \
240  sum += (xx[k] - tmp); \
241  tmp = tmp + sum / n; \
242  } \
243  } \
244  _X_##m [i] = (double)tmp; \
245  }
246 
247 
248 static void
249 cov_complete1(int n, int ncx, double *x, double *xm,
250  int *ind, double *ans, Rboolean *sd_0, Rboolean cor,
251  Rboolean kendall)
252 {
253  COV_init(ncx);
254 
255  if(!kendall) {
256  MEAN(x);/* -> xm[] */
257  n1 = nobs - 1;
258  }
259  for (i = 0 ; i < ncx ; i++) {
260  xx = &x[i * n];
261 
262  if(!kendall) {
263  xxm = xm[i];
264  for (j = 0 ; j <= i ; j++) {
265  yy = &x[j * n];
266  yym = xm[j];
267  sum = 0.;
268  for (k = 0 ; k < n ; k++)
269  if (ind[k] != 0)
270  sum += (xx[k] - xxm) * (yy[k] - yym);
271  ANS(j,i) = ANS(i,j) = (double)(sum / n1);
272  }
273  }
274  else { /* Kendall's tau */
275  for (j = 0 ; j <= i ; j++) {
276  yy = &x[j * n];
277  sum = 0.;
278  for (k = 0 ; k < n ; k++)
279  if (ind[k] != 0)
280  for (n1 = 0 ; n1 < n ; n1++)
281  if (ind[n1] != 0)
282  sum += sign(xx[k] - xx[n1])
283  * sign(yy[k] - yy[n1]);
284  ANS(j,i) = ANS(i,j) = (double)sum;
285  }
286  }
287  }
288 
289  if (cor) {
290  for (i = 0 ; i < ncx ; i++)
291  xm[i] = sqrt(ANS(i,i));
292  for (i = 0 ; i < ncx ; i++) {
293  for (j = 0 ; j < i ; j++) {
294  if (xm[i] == 0 || xm[j] == 0) {
295  *sd_0 = TRUE;
296  ANS(j,i) = ANS(i,j) = NA_REAL;
297  }
298  else {
299  sum = ANS(i,j) / (xm[i] * xm[j]);
300  if(sum > 1.) sum = 1.;
301  ANS(j,i) = ANS(i,j) = (double)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  if(sum > 1.) sum = 1.;
371  ANS(j,i) = ANS(i,j) = (double)sum;
372  }
373  }
374  ANS(i,i) = 1.0;
375  }
376  }
377 } /* cov_na_1() */
378 
379 static void
380 cov_complete2(int n, int ncx, int ncy, double *x, double *y,
381  double *xm, double *ym, int *ind,
382  double *ans, Rboolean *sd_0, Rboolean cor, Rboolean kendall)
383 {
384  COV_init(ncy);
385 
386  if(!kendall) {
387  MEAN(x);/* -> xm[] */
388  MEAN(y);/* -> ym[] */
389  n1 = nobs - 1;
390  }
391  for (i = 0 ; i < ncx ; i++) {
392  xx = &x[i * n];
393  if(!kendall) {
394  xxm = xm[i];
395  for (j = 0 ; j < ncy ; j++) {
396  yy = &y[j * n];
397  yym = ym[j];
398  sum = 0.;
399  for (k = 0 ; k < n ; k++)
400  if (ind[k] != 0)
401  sum += (xx[k] - xxm) * (yy[k] - yym);
402  ANS(i,j) = (double)(sum / n1);
403  }
404  }
405  else { /* Kendall's tau */
406  for (j = 0 ; j < ncy ; j++) {
407  yy = &y[j * n];
408  sum = 0.;
409  for (k = 0 ; k < n ; k++)
410  if (ind[k] != 0)
411  for (n1 = 0 ; n1 < n ; n1++)
412  if (ind[n1] != 0)
413  sum += sign(xx[k] - xx[n1])
414  * sign(yy[k] - yy[n1]);
415  ANS(i,j) = (double)sum;
416  }
417  }
418  }
419 
420  if (cor) {
421 
422 #define COV_SDEV(_X_) \
423  for (i = 0 ; i < nc##_X_ ; i++) { /* Var(X[i]) */ \
424  xx = &_X_[i * n]; \
425  sum = 0.; \
426  if(!kendall) { \
427  xxm = _X_##m [i]; \
428  for (k = 0 ; k < n ; k++) \
429  if (ind[k] != 0) \
430  sum += (xx[k] - xxm) * (xx[k] - xxm); \
431  sum /= n1; \
432  } \
433  else { /* Kendall's tau */ \
434  for (k = 0 ; k < n ; k++) \
435  if (ind[k] != 0) \
436  for (n1 = 0 ; n1 < n ; n1++) \
437  if (ind[n1] != 0 && xx[k] != xx[n1]) \
438  sum ++; /* = sign(. - .)^2 */ \
439  } \
440  _X_##m [i] = (double)SQRTL(sum); \
441  }
442 
443  COV_SDEV(x); /* -> xm[.] */
444  COV_SDEV(y); /* -> ym[.] */
445 
446  for (i = 0 ; i < ncx ; i++)
447  for (j = 0 ; j < ncy ; j++)
448  if (xm[i] == 0. || ym[j] == 0.) {
449  *sd_0 = TRUE;
450  ANS(i,j) = NA_REAL;
451  }
452  else {
453  ANS(i,j) /= (xm[i] * ym[j]);
454  if(ANS(i,j) > 1.) ANS(i,j) = 1.;
455  }
456  }/* cor */
457 
458 }/* cov_complete2 */
459 #undef COV_SDEV
460 
461 static void
462 cov_na_2(int n, int ncx, int ncy, double *x, double *y,
463  double *xm, double *ym, int *has_na_x, int *has_na_y,
464  double *ans, Rboolean *sd_0, Rboolean cor, Rboolean kendall)
465 {
466  COV_ini_na(ncy);
467 
468  if(!kendall) {
469  MEAN_(x, has_na_x);/* -> xm[] */
470  MEAN_(y, has_na_y);/* -> ym[] */
471  n1 = n - 1;
472  }
473  for (i = 0 ; i < ncx ; i++) {
474  if(has_na_x[i]) {
475  for (j = 0 ; j < ncy; j++)
476  ANS(i,j) = NA_REAL;
477  }
478  else {
479  xx = &x[i * n];
480  if(!kendall) {
481  xxm = xm[i];
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  yym = ym[j];
488  sum = 0.;
489  for (k = 0 ; k < n ; k++)
490  sum += (xx[k] - xxm) * (yy[k] - yym);
491  ANS(i,j) = (double)(sum / n1);
492  }
493  }
494  else { /* Kendall's tau */
495  for (j = 0 ; j < ncy ; j++)
496  if(has_na_y[j]) {
497  ANS(i,j) = NA_REAL;
498  } else {
499  yy = &y[j * n];
500  sum = 0.;
501  for (k = 0 ; k < n ; k++)
502  for (n1 = 0 ; n1 < n ; n1++)
503  sum += sign(xx[k] - xx[n1]) * sign(yy[k] - yy[n1]);
504  ANS(i,j) = (double)sum;
505  }
506  }
507  }
508  }
509 
510  if (cor) {
511 
512 #define COV_SDEV(_X_) \
513  for (i = 0 ; i < nc##_X_ ; i++) \
514  if(!has_na_##_X_[i]) { /* Var(X[j]) */ \
515  xx = &_X_[i * n]; \
516  sum = 0.; \
517  if(!kendall) { \
518  xxm = _X_##m [i]; \
519  for (k = 0 ; k < n ; k++) \
520  sum += (xx[k] - xxm) * (xx[k] - xxm); \
521  sum /= n1; \
522  } \
523  else { /* Kendall's tau */ \
524  for (k = 0 ; k < n ; k++) \
525  for (n1 = 0 ; n1 < n ; n1++) \
526  if (xx[k] != xx[n1]) \
527  sum ++; /* = sign(. - .)^2 */ \
528  } \
529  _X_##m [i] = (double) SQRTL(sum); \
530  }
531 
532  COV_SDEV(x); /* -> xm[.] */
533  COV_SDEV(y); /* -> ym[.] */
534 
535  for (i = 0 ; i < ncx ; i++)
536  if(!has_na_x[i]) {
537  for (j = 0 ; j < ncy ; j++)
538  if(!has_na_y[j]) {
539  if (xm[i] == 0. || ym[j] == 0.) {
540  *sd_0 = TRUE;
541  ANS(i,j) = NA_REAL;
542  }
543  else {
544  ANS(i,j) /= (xm[i] * ym[j]);
545  if(ANS(i,j) > 1.) ANS(i,j) = 1.;
546  }
547  }
548  }
549  }/* cor */
550 
551 }/* cov_na_2 */
552 
553 #undef ANS
554 #undef COV_init
555 #undef MEAN
556 #undef MEAN_
557 #undef COV_SDEV
558 
559 /* complete[12]() returns indicator vector ind[] of complete.cases(), or
560  * -------------- if(na_fail) signals error if any NA/NaN is encountered
561  */
562 
563 /* This might look slightly inefficient, but it is designed to
564  * optimise paging in virtual memory systems ...
565  * (or at least that's my story, and I'm sticking to it.)
566 */
567 #define NA_LOOP \
568  for (i = 0 ; i < n ; i++) \
569  if (ISNAN(z[i])) { \
570  if (na_fail) error(_("missing observations in cov/cor"));\
571  else ind[i] = 0; \
572  }
573 
574 #define COMPLETE_1 \
575  double *z; \
576  int i, j; \
577  for (i = 0 ; i < n ; i++) \
578  ind[i] = 1; \
579  for (j = 0 ; j < ncx ; j++) { \
580  z = &x[j * n]; \
581  NA_LOOP \
582  }
583 
584 static void complete1(int n, int ncx, double *x, int *ind, Rboolean na_fail)
585 {
586  COMPLETE_1
587 }
588 
589 static void
590 complete2(int n, int ncx, int ncy, double *x, double *y, int *ind, Rboolean na_fail)
591 {
592  COMPLETE_1
593 
594  for(j = 0 ; j < ncy ; j++) {
595  z = &y[j * n];
596  NA_LOOP
597  }
598 }
599 
600 #define HAS_NA_1(_X_,_HAS_NA_) \
601  for (j = 0 ; j < nc##_X_ ; j++) { \
602  z = &_X_[j * n]; \
603  _HAS_NA_[j] = 0; \
604  for (i = 0 ; i < n ; i++) \
605  if (ISNAN(z[i])) { \
606  _HAS_NA_[j] = 1; break; \
607  } \
608  }
609 
610 
611 static void find_na_1(int n, int ncx, double *x, int *has_na)
612 {
613  double *z;
614  int i, j;
615  HAS_NA_1(x, has_na)
616 }
617 
618 static void
619 find_na_2(int n, int ncx, int ncy, double *x, double *y, int *has_na_x, int *has_na_y)
620 {
621  double *z;
622  int i, j;
623  HAS_NA_1(x, has_na_x)
624  HAS_NA_1(y, has_na_y)
625 }
626 
627 #undef NA_LOOP
628 #undef COMPLETE_1
629 #undef HAS_NA_1
630 
631 /* co[vr](x, y, use =
632  { 1, 2, 3, 4, 5 }
633  "all.obs", "complete.obs", "pairwise.complete", "everything", "na.or.complete"
634  kendall = TRUE/FALSE)
635 */
636 static SEXP corcov(SEXP x, SEXP y, SEXP na_method, SEXP skendall, Rboolean cor)
637 {
638  SEXP ans, xm, ym, ind;
639  Rboolean ansmat, kendall, pair, na_fail, everything, sd_0, empty_err;
640  int i, method, n, ncx, ncy, nprotect = 2;
641 
642  /* Arg.1: x */
643  if(isNull(x)) /* never allowed */
644  error(_("'x' is NULL"));
645 #ifdef _R_in_2017_
646  if(isFactor(x)) error(_("'x' is a factor"));
647 #else
648 # 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."
649  if(isFactor(x)) warning(_(VAR_FACTOR_MSG));
650 #endif
651  /* length check of x -- only if(empty_err) --> below */
652  x = PROTECT(coerceVector(x, REALSXP));
653  if ((ansmat = isMatrix(x))) {
654  n = nrows(x);
655  ncx = ncols(x);
656  }
657  else {
658  n = length(x);
659  ncx = 1;
660  }
661  /* Arg.2: y */
662  if (isNull(y)) {/* y = x : var() */
663  ncy = ncx;
664  } else {
665 #ifdef _R_in_2017_
666  if(isFactor(y)) error(_("'y' is a factor"));
667 #else
669 #endif
671  nprotect++;
672  if (isMatrix(y)) {
673  if (nrows(y) != n)
674  error(_("incompatible dimensions"));
675  ncy = ncols(y);
676  ansmat = TRUE;
677  }
678  else {
679  if (length(y) != n)
680  error(_("incompatible dimensions"));
681  ncy = 1;
682  }
683  }
684  /* Arg.3: method */
685  method = asInteger(na_method);
686 
687  /* Arg.4: kendall */
688  kendall = asLogical(skendall);
689 
690  /* "default: complete" (easier for -Wall) */
691  na_fail = FALSE; everything = FALSE; empty_err = TRUE;
692  pair = FALSE;
693  switch(method) {
694  case 1: /* use all : no NAs */
695  na_fail = TRUE;
696  break;
697  case 2: /* complete */
698  /* did na.omit in R */
699  if (!LENGTH(x)) error(_("no complete element pairs"));
700  break;
701  case 3: /* pairwise.complete */
702  pair = TRUE;
703  break;
704  case 4: /* "everything": NAs are propagated */
705  everything = TRUE;
706  empty_err = FALSE;
707  break;
708  case 5: /* "na.or.complete": NAs are propagated */
709  empty_err = FALSE;
710  break;
711  default:
712  error(_("invalid 'use' (computational method)"));
713  }
714  if (empty_err && !LENGTH(x))
715  error(_("'x' is empty"));
716 
717  if (ansmat) PROTECT(ans = allocMatrix(REALSXP, ncx, ncy));
718  else PROTECT(ans = allocVector(REALSXP, ncx * ncy));
719  sd_0 = FALSE;
720  if (isNull(y)) {
721  if (everything) { /* NA's are propagated */
722  PROTECT(xm = allocVector(REALSXP, ncx));
723  PROTECT(ind = allocVector(LGLSXP, ncx));
724  find_na_1(n, ncx, REAL(x), /* --> has_na[] = */ LOGICAL(ind));
725  cov_na_1 (n, ncx, REAL(x), REAL(xm), LOGICAL(ind), REAL(ans), &sd_0, cor, kendall);
726 
727  UNPROTECT(2);
728  }
729  else if (!pair) { /* all | complete "var" */
730  PROTECT(xm = allocVector(REALSXP, ncx));
731  PROTECT(ind = allocVector(INTSXP, n));
732  complete1(n, ncx, REAL(x), INTEGER(ind), na_fail);
733  cov_complete1(n, ncx, REAL(x), REAL(xm),
734  INTEGER(ind), REAL(ans), &sd_0, cor, kendall);
735  if(empty_err) {
736  Rboolean indany = FALSE;
737  for(i = 0; i < n; i++) {
738  if(INTEGER(ind)[i] == 1) { indany = TRUE; break; }
739  }
740  if(!indany) error(_("no complete element pairs"));
741  }
742  UNPROTECT(2);
743  }
744  else { /* pairwise "var" */
745  cov_pairwise1(n, ncx, REAL(x), REAL(ans), &sd_0, cor, kendall);
746  }
747  }
748  else { /* Co[vr] (x, y) */
749  if (everything) {
750  SEXP has_na_y;
751  PROTECT(xm = allocVector(REALSXP, ncx));
752  PROTECT(ym = allocVector(REALSXP, ncy));
753  PROTECT(ind = allocVector(LGLSXP, ncx));
754  PROTECT(has_na_y = allocVector(LGLSXP, ncy));
755 
756  find_na_2(n, ncx, ncy, REAL(x), REAL(y), INTEGER(ind), INTEGER(has_na_y));
757  cov_na_2 (n, ncx, ncy, REAL(x), REAL(y), REAL(xm), REAL(ym),
758  INTEGER(ind), INTEGER(has_na_y), REAL(ans), &sd_0, cor, kendall);
759  UNPROTECT(4);
760  }
761  else if (!pair) { /* all | complete */
762  PROTECT(xm = allocVector(REALSXP, ncx));
763  PROTECT(ym = allocVector(REALSXP, ncy));
764  PROTECT(ind = allocVector(INTSXP, n));
765  complete2(n, ncx, ncy, REAL(x), REAL(y), INTEGER(ind), na_fail);
766  cov_complete2(n, ncx, ncy, REAL(x), REAL(y), REAL(xm), REAL(ym),
767  INTEGER(ind), REAL(ans), &sd_0, cor, kendall);
768  if(empty_err) {
769  Rboolean indany = FALSE;
770  for(i = 0; i < n; i++) {
771  if(INTEGER(ind)[i] == 1) { indany = TRUE; break; }
772  }
773  if(!indany) error(_("no complete element pairs"));
774  }
775  UNPROTECT(3);
776  }
777  else { /* pairwise */
778  cov_pairwise2(n, ncx, ncy, REAL(x), REAL(y), REAL(ans),
779  &sd_0, cor, kendall);
780  }
781  }
782  if (ansmat) { /* set dimnames() when applicable */
783  if (isNull(y)) {
784  x = getAttrib(x, R_DimNamesSymbol);
785  if (!isNull(x) && !isNull(VECTOR_ELT(x, 1))) {
786  PROTECT(ind = allocVector(VECSXP, 2));
787  SET_VECTOR_ELT(ind, 0, duplicate(VECTOR_ELT(x, 1)));
788  SET_VECTOR_ELT(ind, 1, duplicate(VECTOR_ELT(x, 1)));
789  setAttrib(ans, R_DimNamesSymbol, ind);
790  UNPROTECT(1);
791  }
792  }
793  else {
794  x = getAttrib(x, R_DimNamesSymbol);
796  if ((length(x) >= 2 && !isNull(VECTOR_ELT(x, 1))) ||
797  (length(y) >= 2 && !isNull(VECTOR_ELT(y, 1)))) {
798  PROTECT(ind = allocVector(VECSXP, 2));
799  if (length(x) >= 2 && !isNull(VECTOR_ELT(x, 1)))
800  SET_VECTOR_ELT(ind, 0, duplicate(VECTOR_ELT(x, 1)));
801  if (length(y) >= 2 && !isNull(VECTOR_ELT(y, 1)))
802  SET_VECTOR_ELT(ind, 1, duplicate(VECTOR_ELT(y, 1)));
803  setAttrib(ans, R_DimNamesSymbol, ind);
804  UNPROTECT(1);
805  }
806  }
807  }
808  if(sd_0)/* only in cor() */
809  warning(_("the standard deviation is zero"));
810  UNPROTECT(nprotect);
811  return ans;
812 }
#define VECSXP
Definition: Rinternals.h:124
#define NA_REAL
Definition: Arith.h:65
#define UNPROTECT(n)
Definition: Rinternals.h:653
LibExtern SEXP R_DimNamesSymbol
Definition: Rinternals.h:694
double *() REAL(SEXP x)
Definition: memory.c:3415
Rboolean
Definition: Boolean.h:31
#define REALSXP
Definition: Rinternals.h:117
SEXP SET_VECTOR_ELT(SEXP x, R_xlen_t i, SEXP v)
Definition: memory.c:3451
#define MEAN_(_X_, _HAS_NA_)
Definition: cov.c:226
SEXP coerceVector(SEXP v, SEXPTYPE type)
Definition: coerce.c:1141
#define HAS_NA_1(_X_, _HAS_NA_)
Definition: cov.c:600
void warning(const char *format,...)
Definition: errors.c:262
SEXP() VECTOR_ELT(SEXP x, R_xlen_t i)
Definition: memory.c:3383
j
Definition: qsort-body.c:133
lzma_index ** i
Definition: index.h:625
#define COV_init(_ny_)
Definition: cov.c:190
#define y(k)
#define _(String)
Definition: cov.c:40
INLINE_FUN R_len_t length(SEXP s)
Definition: Rinlinedfuns.h:122
INLINE_FUN Rboolean isMatrix(SEXP s)
Definition: Rinlinedfuns.h:495
static void n
Definition: printvector.c:263
void error(const char *format,...)
Definition: errors.c:758
int *() LOGICAL(SEXP x)
Definition: memory.c:3393
SEXP setAttrib(SEXP vec, SEXP name, SEXP val)
Definition: attrib.c:214
Definition: inflate.h:48
SEXP cor(SEXP x, SEXP y, SEXP na_method, SEXP kendall)
Definition: cov.c:46
int asLogical(SEXP x)
Definition: coerce.c:1602
int nrows(SEXP s)
Definition: util.c:76
int ncols(SEXP s)
Definition: util.c:92
#define COMPLETE_1
Definition: cov.c:574
#define MEAN(_X_)
Definition: cov.c:206
#define COV_SDEV(_X_)
#define isNull
Definition: Rinternals.h:1177
#define ANS(I, J)
Definition: cov.c:64
int asInteger(SEXP x)
Definition: coerce.c:1631
#define COV_ini_na(_ny_)
Definition: cov.c:200
#define TRUE
INTt k
Definition: qsort-body.c:48
#define NA_LOOP
Definition: cov.c:567
#define COV_PAIRWISE_BODY
Definition: cov.c:73
struct SEXPREC * SEXP
Definition: Rinternals.h:438
#define PROTECT(s)
Definition: Rinternals.h:652
#define FALSE
#define LGLSXP
Definition: Rinternals.h:114
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:532
SEXP duplicate(SEXP s)
Definition: duplicate.c:126
INLINE_FUN SEXP allocVector(SEXPTYPE type, R_xlen_t length)
Definition: Rinlinedfuns.h:187
#define INTSXP
Definition: Rinternals.h:116
SEXP allocMatrix(SEXPTYPE mode, int nrow, int ncol)
Definition: array.c:199
#define VAR_FACTOR_MSG
int *() INTEGER(SEXP x)
Definition: memory.c:3401