#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() */