/* $Id: cpl_tools_body.h,v 1.1 2007/04/26 09:17:16 cizzo Exp $
 *
 * This file is part of the ESO Common Pipeline Library
 * Copyright (C) 2001-2004 European Southern Observatory
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 */

/*
 * $Author: cizzo $
 * $Date: 2007/04/26 09:17:16 $
 * $Revision: 1.1 $
 * $Name:  $
 */

/*----------------------------------------------------------------------------*/
/**
  @brief    Sort an array
  @param    pix_arr      array to sort
  @param    n            array size
  @param    reverse      flag indicating whether to sort ascending (zero) or
                         descending (non-zero)
  @param    stable       flag to indicate whether to guarantee stability at
                         the cost of execution time (if cpl_tools_sort_int is
                         O(n log n), this function would still be O(n log n))
  @param    sort_pattern resulting sort pattern
  @return   the #_cpl_error_code_ or CPL_ERROR_NONE

  The heap sort algorithm used here
  - has bad cache performance, but
  - is worst case O(n log n)
  - is stable (meaning that the order of equal elements is conserved 
               indepent on the value of the reverse flag)
  - is in place

  Possible #_cpl_error_code_ set in this function:
  - CPL_ERROR_NULL_INPUT
  - CPL_ERROR_ILLEGAL_INPUT
 */
/*----------------------------------------------------------------------------*/
cpl_error_code
CONCAT2X(cpl_tools_sort_stable_pattern, CPL_TYPE_NAME)(
    CPL_TYPE const *pix_arr,
    int n,
    int reverse,
    int stable,
    int *sort_pattern)
{
    int i;  
    
    /* Check entries */
    cpl_ensure(pix_arr, CPL_ERROR_NULL_INPUT,
               CPL_ERROR_NULL_INPUT) ;
    
    cpl_ensure(sort_pattern, CPL_ERROR_NULL_INPUT,
               CPL_ERROR_NULL_INPUT) ;

    if (n == 0) return CPL_ERROR_NONE;

    for (i = 0; i < n; i++) {
        sort_pattern[i] = i;
    }

    /* 
     * Heap sort
     */
    for (i = n / 2 - 1; i >= 0; i--) {
        int done = 0;
        int root = i;
        int bottom = n - 1;
        
        while ((root*2 + 1 <= bottom) && (!done)) {
            int child = root*2 + 1;

            if (child+1 <= bottom) {
                if ((!reverse && CPL_TOOLS_SORT_LT(
                         pix_arr[sort_pattern[child]],
                         pix_arr[sort_pattern[child + 1]]))
                    ||
                    (reverse && CPL_TOOLS_SORT_LT(
                        pix_arr[sort_pattern[child + 1]],
                        pix_arr[sort_pattern[child]]))
                    ) {
                    child += 1;
                }
            }

            if ((!reverse && CPL_TOOLS_SORT_LT(
                    pix_arr[sort_pattern[root]],
                    pix_arr[sort_pattern[child]])) 
                    ||
                    (reverse && CPL_TOOLS_SORT_LT(
                    pix_arr[sort_pattern[child]],
                    pix_arr[sort_pattern[root]]))
                    ) {
                CPL_INT_SWAP(sort_pattern[root], sort_pattern[child])
                    root = child;
            }
            else {
                done = 1;
            }
        }
    }
    
    for (i = n - 1; i >= 1; i--) {
        int done = 0;
        int root = 0;
        int bottom = i - 1;
        CPL_INT_SWAP(sort_pattern[0], sort_pattern[i])
        
        while ((root*2 + 1 <= bottom) && (!done)) {
            int child = root*2 + 1;

            if (child+1 <= bottom) {
                if ((!reverse && CPL_TOOLS_SORT_LT(
                         pix_arr[sort_pattern[child]],
                         pix_arr[sort_pattern[child + 1]]))
                    ||
                    (reverse && CPL_TOOLS_SORT_LT(
                        pix_arr[sort_pattern[child + 1]],
                        pix_arr[sort_pattern[child]]))
                    ) {
                    child += 1;
                }
            }
            if ((!reverse && CPL_TOOLS_SORT_LT(
                    pix_arr[sort_pattern[root]],
                    pix_arr[sort_pattern[child]])) 
                    ||
                (reverse && CPL_TOOLS_SORT_LT(
                    pix_arr[sort_pattern[child]],
                    pix_arr[sort_pattern[root]]))
                ) {
                CPL_INT_SWAP(sort_pattern[root], sort_pattern[child])
                    root = child;
            }
            else {
                done = 1;
            }
        }
    }

    /* 
     * Enforce stability
     */
    if (stable) {
        for (i = 0; i < n; i++) {
            int j;
            j = i + 1;
            while(j < n &&
                  !CPL_TOOLS_SORT_LT(pix_arr[sort_pattern[i]],
                                     pix_arr[sort_pattern[j]]) &&
                  !CPL_TOOLS_SORT_LT(pix_arr[sort_pattern[j]],
                                     pix_arr[sort_pattern[i]])) {
                j++;
            }
            if (j - i > 1) {
                cpl_tools_sort_int(sort_pattern + i, j - i);
            }
            i = j - 1; 
        }
    }     

    return CPL_ERROR_NONE ;
}