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