"Fossies" - the Fresh Open Source Software Archive

Member "gnuastro-0.8/lib/permutation.c" (3 May 2018, 5423 Bytes) of package /linux/privat/gnuastro-0.8.tar.lz:


As a special service "Fossies" has tried to format the requested source page into HTML format using (guessed) C and C++ source code syntax highlighting (style: standard) with prefixed line numbers and code folding option. Alternatively you can here view or download the uninterpreted source code file. For more information about "permutation.c" see the Fossies "Dox" file reference documentation and the last Fossies "Diffs" side-by-side code changes report: 0.5_vs_0.6.

    1 /*********************************************************************
    2 Permutation -- Work on permutations (arrays of indexs).
    3 This is part of GNU Astronomy Utilities (Gnuastro) package.
    4 
    5 Original author:
    6      Mohammad Akhlaghi <mohammad@akhlaghi.org>
    7 Contributing author(s):
    8 Copyright (C) 2017-2018, Free Software Foundation, Inc.
    9 
   10 Gnuastro is free software: you can redistribute it and/or modify it
   11 under the terms of the GNU General Public License as published by the
   12 Free Software Foundation, either version 3 of the License, or (at your
   13 option) any later version.
   14 
   15 Gnuastro is distributed in the hope that it will be useful, but
   16 WITHOUT ANY WARRANTY; without even the implied warranty of
   17 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
   18 General Public License for more details.
   19 
   20 You should have received a copy of the GNU General Public License
   21 along with Gnuastro. If not, see <http://www.gnu.org/licenses/>.
   22 **********************************************************************/
   23 #include <config.h>
   24 
   25 #include <stdio.h>
   26 #include <string.h>
   27 #include <stdlib.h>
   28 
   29 #include <gnuastro/pointer.h>
   30 #include <gnuastro/permutation.h>
   31 
   32 
   33 
   34 /*********************************************************************/
   35 /***************          Permutation info         *******************/
   36 /*********************************************************************/
   37 void
   38 gal_permutation_check(size_t *permutation, size_t size)
   39 {
   40   size_t i;
   41   for(i=0; i<size; ++i)
   42     printf("after[ %-5zu ]    =   before [ %-5zu ]\n", i, permutation[i]);
   43 }
   44 
   45 
   46 
   47 
   48 
   49 
   50 
   51 
   52 
   53 
   54 
   55 
   56 
   57 
   58 
   59 
   60 
   61 
   62 
   63 
   64 /*********************************************************************/
   65 /***************          Apply permutation        *******************/
   66 /*********************************************************************/
   67 /* Re-order the input dataset based on the given permutation. If
   68    `permutation' is NULL, then the input won't be touched (no re-ordering).
   69 
   70    This is a re-implementation of GSL's `gsl_permute' function (from its
   71    `permutation/permute_source.c'). The reason we didn't use that function
   72    was that it uses system-specific types (like `long' and `int') which are
   73    not easily convertable to Gnuastro's width-based types. There is also a
   74    separate function for each type, heavily using macros to allow a "base"
   75    function to work on all the types. Thus it is hard to
   76    read/understand. Since we use fixed-width types, we can easily use
   77    `memcpy' and have a type-agnostic implementation (only needing the width
   78    of the type).
   79 
   80    As described in GSL's source code and manual, this implementation comes
   81    from Knuth's "Art of computer programmin" book, the "Sorting and
   82    Searching" chapter of Volume 3 (3rd ed). Section 5.2 Exercise 10
   83    (answers), p 617. See there fore more explanations. The algorithm is a
   84    little too abstract, but memory and CPU efficient.
   85 
   86    Definition of permutations:
   87 
   88       permute:    OUT[ i       ]   =   IN[ perm[i] ]     i = 0 .. N-1
   89       inverse:    OUT[ perm[i] ]   =   IN[ i       ]     i = 0 .. N-1
   90 */
   91 void
   92 gal_permutation_apply(gal_data_t *input, size_t *permutation)
   93 {
   94   void *tmp;
   95   size_t i, k, pk, width;
   96   uint8_t *array=input->array;
   97 
   98   /* If permutation is NULL, then it is assumed that the data doesn't need
   99      to be re-ordered. */
  100   if(permutation)
  101     {
  102       /* Necessary initializations. */
  103       width=gal_type_sizeof(input->type);
  104       tmp=gal_pointer_allocate(input->type, 1, 0, __func__, "tmp");
  105 
  106       /* Do the permutation. */
  107       for(i=0;i<input->size;++i)
  108         {
  109           k=permutation[i];
  110 
  111           while(k>i) k=permutation[k];
  112 
  113           if(k>=i)
  114             {
  115               pk = permutation[k];
  116               if( pk != i )
  117                 {
  118                   memcpy(tmp, &array[i*width], width);
  119 
  120                   while(pk!=i)
  121                     {
  122                       memcpy(&array[k*width], &array[pk*width], width);
  123                       k  = pk;
  124                       pk = permutation[k];
  125                     }
  126 
  127                   memcpy(&array[k*width], tmp, width);
  128                 }
  129             }
  130         }
  131 
  132       /* Clean up. */
  133       free(tmp);
  134     }
  135 }
  136 
  137 
  138 
  139 
  140 
  141 /* Apply the inverse of given permutation on the input dataset, see
  142    `gal_permutation_apply_inverse'. */
  143 void
  144 gal_permutation_apply_inverse(gal_data_t *input, size_t *permutation)
  145 {
  146   void *tmp, *ttmp;
  147   size_t i, k, pk, width;
  148   uint8_t *array=input->array;
  149 
  150   if(permutation)
  151     {
  152       /* Initializations */
  153       width=gal_type_sizeof(input->type);
  154       tmp=gal_pointer_allocate(input->type, 1, 0, __func__, "tmp");
  155       ttmp=gal_pointer_allocate(input->type, 1, 0, __func__, "ttmp");
  156 
  157       /* Re-order the values. */
  158       for(i=0;i<input->size;++i)
  159         {
  160           k=permutation[i];
  161 
  162           while(k>i) k=permutation[k];
  163 
  164           if(k>=i)
  165             {
  166               pk = permutation[k];
  167 
  168               if(pk!=i)
  169                 {
  170                   memcpy(tmp, &array[k*width], width);
  171 
  172                   while(pk!=i)
  173                     {
  174                       memcpy(ttmp, &array[pk*width], width);
  175                       memcpy(&array[pk*width], tmp, width);
  176                       memcpy(tmp, ttmp, width);
  177                       k  = pk;
  178                       pk = permutation[k];
  179                     }
  180 
  181                   memcpy(&array[pk*width], tmp, width);
  182                 }
  183             }
  184         }
  185 
  186       /* Clean up. */
  187       free(tmp);
  188       free(ttmp);
  189     }
  190 }