;+
; NAME:
; VALUE_LOCATE
;
; AUTHOR:
; Craig B. Markwardt, NASA/GSFC Code 662, Greenbelt, MD 20770
; craigm@lheamail.gsfc.nasa.gov
;
; PURPOSE:
;
; Locate one or more values in a reference array (IDL LE 5.2 compatibility)
;
; CALLING SEQUENCE:
;
; INDICES = VALUE_LOCATE(REF, VALUES)
;
; DESCRIPTION:
;
; VALUE_LOCATE locates the positions of given values within a
; reference array. The reference array need not be regularly
; spaced. This is useful for various searching, sorting and
; interpolation algorithms.
;
; The reference array should be a monotonically increasing or
; decreasing list of values which partition the real numbers. A
; reference array of NBINS numbers partitions the real number line
; into NBINS+1 regions, like so:
;
;
; REF: X[0] X[1] X[2] X[3] X[NBINS-1]
; <----------|-------------|------|---|----...---|--------------->
; INDICES: -1 0 1 2 3 NBINS-1
;
;
; VALUE_LOCATE returns which partition each of the VALUES falls
; into, according to the figure above. For example, a value between
; X[1] and X[2] would return a value of 1. Values below X[0] return
; -1, and above X[NBINS-1] return NBINS-1. Thus, besides the value
; of -1, the returned INDICES refer to the nearest reference value
; to the left of the requested value.
;
; If the reference array is monotonically decreasing then the
; partitions are numbered starting at -1 from the right instead (and
; the returned INDICES refer to the nearest reference value to the
; *right* of the requested value). If the reference array is
; neither monotonically increasing or decreasing the results of
; VALUE_LOCATE are undefined.
;
; VALUE_LOCATE appears as a built-in funcion in IDL v5.3 and later.
; This version of VALUE_LOCATE should work under IDL v4 and later,
; and is intended to provide a portable solution for users who do
; not have the latest version of IDL. The algrorithm in this file
; is slower but not terribly so, than the built-in version.
;
; Users should be able to place this file in their IDL path safely:
; under IDL 5.3 and later, the built-in function will take
; precedence; under IDL 5.2 and earlier, this function will be used.
;
; INPUTS:
;
; REF - the reference array of monotonically increasing or
; decreasing values.
;
; VALUES - a scalar value or array of values to be located in the
; reference array.
;
;
; KEYWORDS:
;
; L64 - (ignored) for compatibility with built-in version.
;
; NO_CROP - if set, and VALUES is outside of the region between X[0]
; and X[NBINS-1], then the returned indices may be *less
; than* -1 or *greater than* NBINS-1. The user is the
; responsible for cropping these values appropriately.
;
; RETURNS:
;
; An array of indices between -1L and NBINS-1. If VALUES is an
; array then the returned array will have the same dimensions.
;
;
; EXAMPLE:
;
; Cast random values into a histogram with bins from 1-10, 10-100,
; 100-1000, and 1000-10,000.
;
; ;; Make bin edges - this is the ref. array
; xbins = 10D^dindgen(5)
;
; ;; Make some random data that ranges from 1 to 10,000
; x = 10D^(randomu(seed,1000)*4)
;
; ;; Find the bin number of each random value
; ii = value_locate(xbins, x)
;
; ;; Histogram the data
; hh = histogram(ii)
;
;
; SEE ALSO:
;
; VALUE_LOCATE (IDL 5.3 and later), HISTOGRAM, CMHISTOGRAM
;
;
; MODIFICATION HISTORY:
; Written and documented, 21 Jan 2001
; Case of XBINS having only one element, CM, 29 Apr 2001
; Handle case of VALUES exactly hitting REF points, CM, 13 Oct 2001
;
; $Id: value_locate.pro,v 1.5 2001/10/13 17:59:34 craigm Exp $
;
;-
; Copyright (C) 2001, 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.
;-
function value_locate, xbins, x, l64=l64, no_crop=nocrop, _EXTRA=extra
on_error, 2
nbins = n_elements(xbins)
sz = size(xbins)
;; Error checking
if nbins EQ 0 then message, 'ERROR: XBINS must have at least one element'
if nbins EQ 1 then return, (x GE xbins(0)) - 1L
;; The values are computed by spline interpolation. Here is the "y"
;; value of the spline, which is just the bin position.
tp = sz(sz(0)+1)
if tp EQ 1 OR tp EQ 2 OR tp EQ 12 then begin
yy = findgen(nbins) - 0.5
eps = (machar()).eps
endif else begin
yy = dindgen(nbins) - 0.5D
eps = (machar(/double)).eps
endelse
;; Check if we are reversing.
if xbins(nbins-1) GT xbins(0) then rev = 0 else rev = 1
;; Compute the spline interpolation. Note here that we set the 2nd
;; derivative value to zero since the derivative computed by
;; SPL_INIT seems to be royally screwed up. Also we do separate
;; computations for the increasing and non-increasing cases, since
;; SPL_INTERP seems to choke on the later.
if rev EQ 0 then $
ii = round(spl_interp(xbins, yy, yy*0, x) + eps) $
else $
ii = round(spl_interp(reverse(xbins), yy, yy*0, x) + eps)
;; Crop the end values appropriately
if NOT keyword_set(nocrop) then begin
ii = ii > (-1L) < (nbins-1)
endif
;; Reverse the array
if rev EQ 0 then $
ret = temporary(ii) $
else $
ret = (nbins-2)-temporary(ii)
;; Reform the array to the correct dimensions (ie, add trailing
;; dimensions)
sz = size(x)
if sz(0) GT 0 then $
ret = reform(ret, sz(1:sz(0)), /overwrite)
return, ret
end