htswanalysis/MACS/lib/gsl/gsl-1.11/permutation/permute_source.c

   1 /* permutation/permute_source.c
   2  *
   3  * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007 Brian Gough
   4  *
   5  * This program is free software; you can redistribute it and/or modify
   6  * it under the terms of the GNU General Public License as published by
   7  * the Free Software Foundation; either version 3 of the License, or (at
   8  * your option) any later version.
   9  *
  10  * This program is distributed in the hope that it will be useful, but
  11  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13  * General Public License for more details.
  14  *
  15  * You should have received a copy of the GNU General Public License
  16  * along with this program; if not, write to the Free Software
  17  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
  18  */
  19
  20 /* In-place Permutations
  21
  22    permute:    OUT[i]       = IN[perm[i]]     i = 0 .. N-1
  23    invpermute: OUT[perm[i]] = IN[i]           i = 0 .. N-1
  24
  25    PERM is an index map, i.e. a vector which contains a permutation of
  26    the integers 0 .. N-1.
  27
  28    From Knuth "Sorting and Searching", Volume 3 (3rd ed), Section 5.2
  29    Exercise 10 (answers), p 617
  30
  31    FIXME: these have not been fully tested.
  32 */
  33
  34 int
  35 TYPE (gsl_permute) (const size_t * p, ATOMIC * data, const size_t stride, const size_t n)
  36 {
  37   size_t i, k, pk;
  38
  39   for (i = 0; i < n; i++)
  40     {
  41       k = p[i];
  42
  43       while (k > i)
  44         k = p[k];
  45
  46       if (k < i)
  47         continue ;
  48
  49       /* Now have k == i, i.e the least in its cycle */
  50
  51       pk = p[k];
  52
  53       if (pk == i)
  54         continue ;
  55
  56       /* shuffle the elements of the cycle */
  57
  58       {
  59         unsigned int a;
  60
  61         ATOMIC t[MULTIPLICITY];
  62
  63         for (a = 0; a < MULTIPLICITY; a++)
  64           t[a] = data[i*stride*MULTIPLICITY + a];
  65
  66         while (pk != i)
  67           {
  68             for (a = 0; a < MULTIPLICITY; a++)
  69               {
  70                 ATOMIC r1 = data[pk*stride*MULTIPLICITY + a];
  71                 data[k*stride*MULTIPLICITY + a] = r1;
  72               }
  73             k = pk;
  74             pk = p[k];
  75           };
  76
  77         for (a = 0; a < MULTIPLICITY; a++)
  78           data[k*stride*MULTIPLICITY + a] = t[a];
  79       }
  80     }
  81
  82   return GSL_SUCCESS;
  83 }
  84
  85 int
  86 FUNCTION (gsl_permute,inverse) (const size_t * p, ATOMIC * data, const size_t stride, const size_t n)
  87 {
  88   size_t i, k, pk;
  89
  90   for (i = 0; i < n; i++)
  91     {
  92       k = p[i];
  93
  94       while (k > i)
  95         k = p[k];
  96
  97       if (k < i)
  98         continue ;
  99
 100       /* Now have k == i, i.e the least in its cycle */
 101
 102       pk = p[k];
 103
 104       if (pk == i)
 105         continue ;
 106
 107       /* shuffle the elements of the cycle in the inverse direction */
 108
 109       {
 110         unsigned int a;
 111
 112         ATOMIC t[MULTIPLICITY];
 113
 114         for (a = 0; a < MULTIPLICITY; a++)
 115           t[a] = data[k*stride*MULTIPLICITY+a];
 116
 117         while (pk != i)
 118           {
 119             for (a = 0; a < MULTIPLICITY; a++)
 120               {
 121                 ATOMIC r1 = data[pk*stride*MULTIPLICITY + a];
 122                 data[pk*stride*MULTIPLICITY + a] = t[a];
 123                 t[a] = r1;
 124               }
 125
 126             k = pk;
 127             pk = p[k];
 128           };
 129
 130         for (a = 0; a < MULTIPLICITY; a++)
 131           data[pk*stride*MULTIPLICITY+a] = t[a];
 132       }
 133     }
 134
 135   return GSL_SUCCESS;
 136 }
 137
 138
 139 int
 140 TYPE (gsl_permute_vector) (const gsl_permutation * p, TYPE (gsl_vector) * v)
 141 {
 142   if (v->size != p->size)
 143     {
 144       GSL_ERROR ("vector and permutation must be the same length", GSL_EBADLEN);
 145     }
 146
 147   TYPE (gsl_permute) (p->data, v->data, v->stride, v->size) ;
 148
 149   return GSL_SUCCESS;
 150 }
 151
 152 int
 153 FUNCTION (gsl_permute_vector,inverse) (const gsl_permutation * p, TYPE (gsl_vector) * v)
 154 {
 155   if (v->size != p->size)
 156     {
 157       GSL_ERROR ("vector and permutation must be the same length", GSL_EBADLEN);
 158     }
 159
 160   FUNCTION (gsl_permute,inverse) (p->data, v->data, v->stride, v->size) ;
 161
 162   return GSL_SUCCESS;
 163 }