;+
; NAME:
;       MAX_CUMULATIVE
;
; PURPOSE:
;
;       Efficiently perform a cumulative MAX over some dimension of an
;       array.
;
; CALLING SEQUENCE:
;
;       max_cumulative, array [, dimension]
;
; INPUTS:
;
;       array: An array of at least 2 dimensions upon which to
;          operate.  Will be modified on output.
;
;       dimension: The dimension along which to perform a cumulative
;          maximum, starting at 1 (1:rows, 2:columns, ...).  Defaults
;          to the final dimension.
;
; OUTPUTS:
;
;       array: the array is modified on output to contain the
;          cumulative maximum along the specified dimension.
;
; EXAMPLE:
;       
;       a=randomu(sd,5,4,3,2)
;       max_cumulative,a,2
;
; SEE ALSO:
;
;       SORT_ND, HIST_ND
;
; MODIFICATION HISTORY:
;
;       Fri Mar 9 11:51:59 2012, J.D. Smith: Removed index array for
;          further speedup (suggestion of H. Stege).
;       Thu Mar 8 12:56:58 2012, J.D. Smith: Simplified
;          pre-transposition and indexing
;       Mon Feb 27 18:41:52 2012, J.D. Smith: Written
;		
;
;-
;##############################################################################
;
; LICENSE
;
;  Copyright (C) 2012 J.D. Smith
;
;  This file 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, or (at your
;  option) any later version.
;
;  This file 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 file; see the file COPYING.  If not, write to the
;  Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
;  Boston, MA 02110-1301, USA.
;
;##############################################################################

pro max_cumulative,a,dim
  s=size(a,/DIMENSIONS)
  ns=n_elements(s) 
  if n_elements(dim) eq 0 then d=ns-1 else d=dim-1
  
  t=where(indgen(ns) ne d)
  if d lt ns-1 then a=transpose(a,[t,d])
  
  off=product(s[t],/PRESERVE_TYPE)
  for i=1,s[d]-1 do a[i*off]=a[(i-1)*off:i*off-1]>a[i*off:(i+1)*off-1]
  
  if d lt ns-1 then a=transpose(a,sort([t,d]))
end