htswanalysis/MACS/lib/gsl/gsl-1.11/sort/subset_source.c

   1 /* sort/subset_source.c
   2  *
   3  * Copyright (C) 1999,2000,2001  Thomas Walter, Brian Gough
   4  *
   5  * This is free software; you can redistribute it and/or modify it
   6  * under the terms of the GNU General Public License as published by the
   7  * Free Software Foundation; either version 3, or (at your option) any
   8  * later version.
   9  *
  10  * This source is distributed in the hope that it will be useful, but WITHOUT
  11  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  12  * FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  13  * for more details.
  14  */
  15
  16 /* find the k-th smallest elements of the vector data, in ascending order */
  17
  18 int
  19 FUNCTION (gsl_sort, smallest) (BASE * dest, const size_t k,
  20                                const BASE * src, const size_t stride,
  21                                const size_t n)
  22 {
  23   size_t i, j;
  24   BASE xbound;
  25
  26   if (k > n)
  27     {
  28       GSL_ERROR ("subset length k exceeds vector length n", GSL_EINVAL);
  29     }
  30
  31   if (k == 0 || n == 0)
  32     {
  33       return GSL_SUCCESS;
  34     }
  35
  36   /* take the first element */
  37
  38   j = 1;
  39   xbound = src[0 * stride];
  40   dest[0] = xbound;
  41
  42   /* examine the remaining elements */
  43
  44   for (i = 1; i < n; i++)
  45     {
  46       size_t i1;
  47
  48       BASE xi = src[i * stride];
  49
  50       if (j < k)
  51         {
  52           j++;
  53         }
  54       else if (xi >= xbound)
  55         {
  56           continue;
  57         }
  58
  59       for (i1 = j - 1; i1 > 0 ; i1--)
  60         {
  61           if (xi > dest[i1 - 1])
  62             break;
  63
  64           dest[i1] = dest[i1 - 1];
  65         }
  66
  67       dest[i1] = xi;
  68
  69       xbound = dest[j-1];
  70     }
  71
  72   return GSL_SUCCESS;
  73 }
  74
  75
  76 int
  77 FUNCTION (gsl_sort_vector,smallest) (BASE * dest, const size_t k,
  78                                      const TYPE (gsl_vector) * v)
  79 {
  80   return FUNCTION (gsl_sort, smallest) (dest, k, v->data, v->stride, v->size);
  81 }
  82
  83 int
  84 FUNCTION (gsl_sort, largest) (BASE * dest, const size_t k,
  85                               const BASE * src, const size_t stride,
  86                               const size_t n)
  87 {
  88   size_t i, j;
  89   BASE xbound;
  90
  91   if (k > n)
  92     {
  93       GSL_ERROR ("subset length k exceeds vector length n", GSL_EINVAL);
  94     }
  95
  96   if (k == 0 || n == 0)
  97     {
  98       return GSL_SUCCESS;
  99     }
 100
 101   /* take the first element */
 102
 103   j = 1;
 104   xbound = src[0 * stride];
 105   dest[0] = xbound;
 106
 107   /* examine the remaining elements */
 108
 109   for (i = 1; i < n; i++)
 110     {
 111       size_t i1;
 112
 113       BASE xi = src[i * stride];
 114
 115       if (j < k)
 116         {
 117           j++;
 118         }
 119       else if (xi <= xbound)
 120         {
 121           continue;
 122         }
 123
 124       for (i1 = j - 1; i1 > 0 ; i1--)
 125         {
 126           if (xi < dest[i1 - 1])
 127             break;
 128
 129           dest[i1] = dest[i1 - 1];
 130         }
 131
 132       dest[i1] = xi;
 133
 134       xbound = dest[j-1];
 135     }
 136
 137   return GSL_SUCCESS;
 138 }
 139
 140
 141 int
 142 FUNCTION (gsl_sort_vector,largest) (BASE * dest, const size_t k,
 143                                     const TYPE (gsl_vector) * v)
 144 {
 145   return FUNCTION (gsl_sort, largest) (dest, k, v->data, v->stride, v->size);
 146 }