1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
|
#include <config.h>
#include <math.h>
#include <gsl/gsl_statistics.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_sort_vector.h>
static int compute_rank(gsl_vector *v);
#define BASE_LONG_DOUBLE
#include "templates_on.h"
#include "covariance_source.c"
#include "templates_off.h"
#undef BASE_LONG_DOUBLE
#define BASE_DOUBLE
#include "templates_on.h"
#include "covariance_source.c"
#include "templates_off.h"
#undef BASE_DOUBLE
#define BASE_FLOAT
#include "templates_on.h"
#include "covariance_source.c"
#include "templates_off.h"
#undef BASE_FLOAT
#define BASE_ULONG
#include "templates_on.h"
#include "covariance_source.c"
#include "templates_off.h"
#undef BASE_ULONG
#define BASE_LONG
#include "templates_on.h"
#include "covariance_source.c"
#include "templates_off.h"
#undef BASE_LONG
#define BASE_UINT
#include "templates_on.h"
#include "covariance_source.c"
#include "templates_off.h"
#undef BASE_UINT
#define BASE_INT
#include "templates_on.h"
#include "covariance_source.c"
#include "templates_off.h"
#undef BASE_INT
#define BASE_USHORT
#include "templates_on.h"
#include "covariance_source.c"
#include "templates_off.h"
#undef BASE_USHORT
#define BASE_SHORT
#include "templates_on.h"
#include "covariance_source.c"
#include "templates_off.h"
#undef BASE_SHORT
#define BASE_UCHAR
#include "templates_on.h"
#include "covariance_source.c"
#include "templates_off.h"
#undef BASE_UCHAR
#define BASE_CHAR
#include "templates_on.h"
#include "covariance_source.c"
#include "templates_off.h"
#undef BASE_CHAR
/*
compute_rank()
Compute rank of a sorted vector
Inputs: v - sorted data vector on input; rank vector on output
Notes: ranks are always computed in double precision
*/
static int
compute_rank(gsl_vector *v)
{
const size_t n = v->size;
size_t i = 0;
while (i < n - 1)
{
double vi = gsl_vector_get(v, i);
if (vi == gsl_vector_get(v, i + 1))
{
size_t j = i + 2;
size_t k;
double rank = 0.0;
/* we have detected a tie, find number of equal elements */
while (j < n && vi == gsl_vector_get(v, j))
++j;
/* compute rank */
for (k = i; k < j; ++k)
rank += k + 1.0;
/* divide by number of ties */
rank /= (double) (j - i);
for (k = i; k < j; ++k)
gsl_vector_set(v, k, rank);
i = j;
}
else
{
/* no tie - set rank to natural ordered position */
gsl_vector_set(v, i, i + 1.0);
++i;
}
}
if (i == n - 1)
gsl_vector_set(v, n - 1, (double) n);
return GSL_SUCCESS;
} /* compute_rank() */
|
