multisort.pro

;+
; NAME:
;   MULTISORT
;
; AUTHOR:
;   Craig B. Markwardt, NASA/GSFC Code 662, Greenbelt, MD 20770
;   craigm@lheamail.gsfc.nasa.gov
;
; PURPOSE:
;   Perform a sorting operation with multiple sort keys
;
; CALLING SEQUENCE:
;   INDICES = MULTISORT(KEY1, KEY2, ..., [/L64, ], [ORDER=order])
;
; DESCRIPTION:
;
;   The function MULTISORT performs a sorting operation with multiple
;   sort keys.  Unlike the IDL built-in SORT() function, which can
;   only sort a single key, MULTISORT can accept multiple keys.  In
;   cases where the primary key is equal, the sort order is based on
;   any secondary keys provided.  The return value is an array of
;   indices which will place the key arrays into sorted order.
;
;   MULTISORT works by building an internal sort key string which can
;   be sorted in a single pass.  Because MULTISORT is not a built-in
;   function, and because it must build these auxiliary strings, it
;   cannot be as fast or memory-efficient as the built-in function.
;   Users will need several times more memory than the memory used
;   to store just the input keys.
;
;   MULTISORT() allows the user to choose the sort order for each key
;   separately.  The ORDER keyword is an N-vector, one order for each
;   input key.  ORDER[i] is +1 to sort KEYi ascending, and ORDER[i] is
;   -1 to sort KEYi descending.
;
; INPUTS:
;
;   KEY1, KEY2, ... KEY9 - input sort keys.  Any integer, floating
;                          point or string value is allowed.  The
;                          number of values must be the same for each
;                          key.
;
;
;
; KEYWORDS:
;
;   ORDER - an N-vector, giving the sort order for each key (see
;           documentation above).
;           Default: REPLICATE(+1,N_PARAMS())  (all keys ascending)
;
;   L64 - if set, then return a LONG64 index instead of LONG.
;
; RETURNS:
;
;   An array of indices which will place the keys into sorted order.
;   I.e., KEYS1[INDICES],  KEYS2[INDICES] ... will be in order.
;
; SEE ALSO:
;
;   SORT
;
; MODIFICATION HISTORY:
;   Written, CM, Jun 2007
;   Document the encoding format, and make some floating point
;     operations more efficient, CM, Jan 2008
;   Correct several bugs in the handling of floating point numbers
;     in the range -1.0 to 1.0, made more efficient, (thanks to Eric
;     Jensen);  I also saved some test cases; CM, Jul 2008
;
;
;  $Id: multisort.pro,v 1.3 2008/07/08 20:21:18 craigm Exp $
;
;-
; Copyright (C) 2007, 2008, Craig Markwardt
; This software is provided as is without any warranty whatsoever.
; Permission to use, copy, modify, and distribute modified or
; unmodified copies is granted, provided this copyright and disclaimer
; are included unchanged.
;-
; ========================================================
; Utility function to transform an integer into a string key
;  Integers are coded in hexidecimal, with a prefix of
;  'A' for negative and 'B' for positive.
function multisort_intkey, x, len, unsigned=u, order=order
  COMPILE_OPT strictarr
  n = n_elements(x)
  if order LT 0 then x1 = NOT temporary(x) $ ;; Reverse order
  else x1 = temporary(x)

  ;; v = [-2L^31L,-32768,-2000,-20,-1,0,1,20,2000,32767,(-2L^31L)-1]

  slen = strtrim(len,2)
  fmt = '(Z'+slen+'.'+slen+')'  ;; (In.n)  - zero-padded
  if NOT keyword_set(u) then begin
    prestr = strarr(n)+'B' ;; Prefix to indicate positive values
    wh = where(x1 LT 0, ct)
    if ct GT 0 then prestr[wh] = 'A'

    bmask = ishft((x1[0] AND 0b) + 1b,4*len)-1
    if bmask LE 1 then bmask = NOT (x1[0] AND 0b)
    x1 = x1 AND bmask

    return, temporary(prestr)+string(x1,format=fmt)
  endif else begin
    return,                   string(x1,format=fmt)
  endelse

end

; ========================================================
; Utility function to transform a float into a string key
; Floating point numbers are converted to IEEE format,
; and then examined as integers, allowing MULTISORT_INTKEY
; to be used.
;
function multisort_fltkey, x1, type, order=order
  COMPILE_OPT strictarr
  n = n_elements(x1)

  if type EQ 4 then begin
    ;; Floating point data (4 bytes)

    ;; Test case
    ;; v = [-!values.f_infinity,-2000.,-20,-1.5,-1.0,-0.5,0,0.5,1.0,1.5,20,2000,!values.f_infinity,!values.f_nan]
    byteorder, x1, /ftoxdr
    x1 = long(temporary(x1),0,n)
    byteorder, x1, /ntohl
    wh = where(x1 LT 0, ct)
    if ct GT 0 then x1[wh] = x1[wh] XOR '7fffffff'xl

    return, multisort_intkey(x1,8,order=order)

  endif else begin
    ;; Double precision data

    ;; Test case
    ;; v = [-!values.d_infinity,-2000d,-20,-1.5,-1.0,-0.5,0,0.5,1.0,1.5,20,2000,!values.d_infinity,!values.d_nan]
    byteorder, x1, /dtoxdr
    x1 = long64(temporary(x1),0,n)
    byteorder, x1, /l64swap, /swap_if_little
    wh = where(x1 LT 0, ct)
    if ct GT 0 then x1[wh] = x1[wh] XOR '7fffffffffffffff'xll

    return, multisort_intkey(x1,16,order=order)

  endelse

end

; ========================================================
; Utility function to transform a string into a string key
function multisort_strkey, x, order=order
  COMPILE_OPT strictarr
  len = strlen(x)
  maxlen = max(len, min=minlen)
  if maxlen GT minlen then begin
    ;; Pad out to the maximum string length (i.e. left-align the strings)
    pad = string(bytarr(maxlen-minlen)+32b)
    key = strmid(x+pad,0,maxlen)
  endif else begin
    key = x
  endelse

  ;; Reverse order if requested
  if order LT 0 then begin
    key = string( (255b - byte(temporary(key))) > 1b )
  endif

  return, key
end


; ========================================================
; MAIN ROUTINE
; ========================================================
function multisort, x0, x1, x2, x3, x4, x5, x6, x7, x8, x9, L64=L64, $
  keys=keys0, order=order0, no_builtin=nobuiltin

  COMPILE_OPT strictarr

  nkeys = n_params()
  if nkeys EQ 0 then begin
    message, 'USAGE: INDICES = MULTISORT(KEY1[,KEY2,KEY3,...])', /info
    return, -1L
  endif

  order = intarr(nkeys) + 1
  if n_elements(order0) GT 0 then order[0] = round(order0)

  ;; Special case: only one term, no need to do complicate sort key
  ;; manipulations.
  if nkeys EQ 1 AND order[0] EQ +1 AND NOT keyword_set(nobuiltin) then begin
    return, sort(x0, L64=L64)
  endif

  ;; Master key
  mkey = ''
  for i = 0, nkeys-1 do begin
    xi = 0 & dummy = temporary(xi)
    case i of
      0: xi = x0
      1: xi = x1
      2: xi = x2
      3: xi = x3
      4: xi = x4
      5: xi = x5
      6: xi = x6
      7: xi = x7
      8: xi = x8
      9: xi = x9
    endcase
    if n_elements(xi) EQ 0 then begin
      message, string(i,format='("ERROR: no data was in parameter X",I0)')
      return, -1L
    endif

    sz = size(xi)
    tp = sz[sz[0]+1]
    o = order[i]
    case tp of
      1:  mkey = temporary(mkey) + multisort_intkey(temporary(xi),2,/u,o=o)  ;; BYTE
      2:  mkey = temporary(mkey) + multisort_intkey(temporary(xi),4,o=o)     ;; INT
      3:  mkey = temporary(mkey) + multisort_intkey(temporary(xi),8,o=o)     ;; LONG
      4:  mkey = temporary(mkey) + multisort_fltkey(temporary(xi),4,o=o)     ;; FLOAT
      5:  mkey = temporary(mkey) + multisort_fltkey(temporary(xi),5,o=o)     ;; DOUBLE
      7:  mkey = temporary(mkey) + multisort_strkey(temporary(xi),o=o)       ;; STRING
      12: mkey = temporary(mkey) + multisort_intkey(temporary(xi),4,/u,o=o)  ;; UINT
      13: mkey = temporary(mkey) + multisort_intkey(temporary(xi),8,/u,o=o)  ;; ULONG
      14: mkey = temporary(mkey) + multisort_intkey(temporary(xi),16,o=o)    ;; LONG64
      15: mkey = temporary(mkey) + multisort_intkey(temporary(xi),16,/u,o=o) ;; ULONG64
      else: begin
        message, string(tp, i, $
          format='("ERROR: data type ",I0," for parameter X,",I0," is not sortable")')
        return, -1L
      end
    endcase

    xi = 0
  endfor

  return, sort(mkey, L64=L64)
end