;
; Test the following problem:
;
;   Given a list of indices with repeats and a data vector of the
;   same length, calculate a result vector of length max(inds)+1
;   such that:
;
;     vec[i]=total(data[where(inds eq i)])
;
;   (or zero if no ind at that point).
;
iter=10
n_ind=50000L
repeat_fac=10L                  ;how many times on average indices are repeated
sparse_fac=1L                   ;how sparsely separated indices are
inds=long(randomu(100,n_ind)*n_ind/repeat_fac*sparse_fac)
;inds=replicate(1L,n_ind)        ;pathological case
data=findgen(n_ind)/5

;Simple method, using where on each index
t=systime(1)
for i=1,iter do begin 
   mx=max(inds)
   vec1=fltarr(mx+1) 
   for j=0L,mx do begin 
      wh=where(inds eq j,cnt)
      if cnt eq 0 then continue
      vec1[j]=total(data[wh])
   endfor 
endfor 
print,FORMAT='("WHERE loop:",T50,F8.4)',(systime(1)-t)/iter

;; Simple method, loop and accumulate
t=systime(1)
for i=1,iter do begin 
   mx=max(inds)
   vec2=fltarr(mx+1) 
   for j=0L,n_elements(data)-1L do $
      vec2[inds[j]]=vec2[inds[j]]+data[j]
endfor 
print,FORMAT='("Literal Accumulate Loop:",T50,F8.4)',(systime(1)-t)/iter

;; Smarter method, looping over reverse indices
t=systime(1)
for i=1,iter do begin 
   mx=max(inds)
   vec3=fltarr(mx+1) 
   h=histogram(inds,reverse_indices=ri,OMIN=om)
   for j=0L,n_elements(h)-1 do if ri[j+1] gt ri[j] then $
      vec3[j+om]=total(data[ri[ri[j]:ri[j+1]-1]])
endfor 
print,FORMAT='("Reverse Indices Loop:",T50,F8.4)',(systime(1)-t)/iter

;; Loop free method using sparse arrays
t=systime(1)
for i=1,iter do begin 
   mx=max(inds)
   h1=histogram(inds,OMIN=om,REVERSE_INDICES=ri1)
   col=ri1[n_elements(h1)+1:*]
   h2=histogram(total(h1,/cumulative)-1,MIN=0,reverse_indices=ri2)
   row=ri2[0:n_elements(h2)-1]-ri2[0]+om ; chunk indices = row number
   sparse_array=sprsin(col,row,replicate(1.,n_ind),(mx+1)>n_ind)
   if mx ge n_ind then $
      vec4=sprsax(sparse_array,[data,replicate(0.,mx+1-n_ind)]) $
   else vec4=(sprsax(sparse_array,data))[0:mx]
endfor 
print,FORMAT='("Loop-Free with Sparse Arrays:",T50,F8.4)',(systime(1)-t)/iter

;; Wayne Landsman's "loop over bin depth" method from FDDRIZZLE
t=systime(1)
for i=1,iter do begin 
   mx=max(inds)
   vec5=fltarr(mx+1)
   h=histogram(inds,REVERSE_INDICES=ri,omin=om)
   gmax = max(h)                ;Highest number of duplicate indicies
   for j=1L,gmax do begin
      g = where(h GE j, Ng)
      if Ng GT 0 then  vec5[om+g] = vec5[om+g] + data[ri[ ri[g]+j-1]]
   endfor
endfor 
print,FORMAT='("FDDRIZZLE Loop:",T50,F8.4)',(systime(1)-t)/iter
   
;; Modified dual histogram "loop over bin depth" method
t=systime(1)
for i=1,iter do begin 
   mx=max(inds)
   vec6=fltarr(mx+1)
   h1=histogram(inds,reverse_indices=ri1,OMIN=om)
   h2=histogram(h1,reverse_indices=ri2,MIN=1)
   ;; easy case - single values w/o duplication
   if ri2[1] gt ri2[0] then begin 
      vec_inds=ri2[ri2[0]:ri2[1]-1] 
      vec6[om+vec_inds]=data[ri1[ri1[vec_inds]]]
   endif
   for j=(min(h1)-1)>1,n_elements(h2)-1 do begin 
      if ri2[j+1] eq ri2[j] then continue ;none with that many duplicates
      vec_inds=ri2[ri2[j]:ri2[j+1]-1] ;indices into h1
      vinds=om+vec_inds
      vec_inds=rebin(ri1[vec_inds],h2[j],j+1,/SAMPLE)+ $
               rebin(transpose(lindgen(j+1)),h2[j],j+1,/SAMPLE)
      vec6[vinds]+=total(data[ri1[vec_inds]],2)
   endfor 
endfor 
print,FORMAT='("Dual Histogram Loop:",T50,F8.4)',(systime(1)-t)/iter

;; Craig Markwardt's FDDRIZZLE-Like thinned WHERE histogram loop.
t=systime(1)
for i=1,iter do begin 
   mx=max(inds)
   vec7=fltarr(mx+1)
   h = histogram(inds,OMIN=om,REVERSE_INDICES=ri)
   wh = where(h GT 0)
   mx = max(h[wh],min=mn)
   for j=1L,mx do begin
      wh=wh[where(h[wh] GE j)]  ; Get IND cells with GE i entries
      vec7[om+wh] += data[ri[ri[wh]+j-1]] ; Add into the total
   endfor
endfor 
print,FORMAT='("Thinned WHERE Histogram Loop:",T50,F8.4)',(systime(1)-t)/iter

;; Sort cumulative total, from Jaco van Gorkom
t=systime(1)
for i=1,iter do begin 
   h = histogram(inds, REVERSE_INDICES=ri)
   nh = n_elements(h)
   tdata = [0.,total(data[ri[nh+1:*]],/CUMULATIVE,/DOUBLE)]
   vec8 = tdata[ri[1:nh]-nh-1]-tdata[ri[0:nh-1]-nh-1]
endfor 
print,FORMAT='("Sort cumulative total:",T50,F8.4)',(systime(1)-t)/iter


END

;; DLM
; t=systime(1)
; for i=1,iter do begin 
;    mx=max(inds)
;    vec9=fltarr(mx+1)
;    chunk_decimate,vec9,inds,data
; endfor 
; print,FORMAT='("Literal Accumulate: Compiled DLM :",T50,F8.4)', $
;       (systime(1)-t)/iter