;
;Copyright 1996-2013 United States Government as represented by the
;Administrator of the National Aeronautics and Space Administration.
;All Rights Reserved.
; This software may be used, copied, or redistributed as long as it is not
; sold and this copyright notice is reproduced on each copy made.  This
; routine is provided as is without any express or implied warranties
; whatsoever.
;
pro cmb_setemptybin2fill,b,b_nbin,fillval
ii=where( b_nbin eq 0) 
if ii[0] ne -1 then b[ii] = fillval
print,'in cmb_setemptybin2fill fillval=',fillval
end


pro cmb_fill_array,b,xb,b_nbin
; fills empty bins by linear interpolating across them.
; input/output
; b[nc,nt]  a two dimensional array of binned data, trailing dimesnion is time, if scaler b[1,nt]
; xb[nt] is center time of the binned data
; exptrapoled bin values set to value of closest valid bin value
; b_nbin[nc,nt] is the number of data points used in each bin to binning.
ii=where(b_nbin eq 0)
if ii[0] eq -1 then return
si0 = size(b,/structure)
ndims = si0.n_dimensions
;help,si0,/str
if ndims eq 1 then begin
        iz = where(b_nbin eq 0) ;find bins with no data
	iv = where(b_nbin ne 0)
	if iz[0] eq -1 or n_elements(iv) lt 2 then return
	y = b
	yz = interpol(b[iv],xb[iv],xb[iz])
	y[iz] =yz
	b[*]=y
	xr = [min(xb[iv],i0),max(xb[iv],i1)]
	i=where(xb lt xr[0]) & if i[0] ne -1 then b[i]= b[iv[i0]]
	i=where(xb gt xr[1]) & if i[0] ne -1 then b[i]= b[iv[i1]] 
        return
endif else begin
    for ic = 0l,si0.dimensions[0]-1 do begin   
	iz = where(b_nbin[ic,*] eq 0) ;find bins with no data
	iv = where(b_nbin[ic,*] ne 0)
	if iz[0] eq -1 or n_elements(iv) lt 2 then goto, next
	y = b[ic,*]
	yz = interpol(b[ic,iv],xb[iv],xb[iz])
	y[iz] =yz
	b[ic,*]=y
	xr = [min(xb[iv],i0),max(xb[iv],i1)]
	i=where(xb lt xr[0]) & if i[0] ne -1 then b[ic,i]= b[ic,iv[i0]]
	i=where(xb gt xr[1]) & if i[0] ne -1 then b[ic,i]= b[ic,iv[i1]] 
	next:
    endfor
endelse
end

pro cmb_timebin_array,tin00,serin0,t0,tf,dtbin $
    ,tout,serout,fillval, seroutuncertainty, missing_variance_flag $
    ,serout_flag=serout_flag, fill_empty_bins=fill_empty_bins $
    ,autobad=autobad,sigmul=sigmul $
   ,VALIDMIN=VALIDMIN, VALIDMAX=VALIDMAX
;input 
; t0-start time of binning
; tf-end time of binning
; dtbin time increment
; note: t0, tf, dtbin, tin0 must be in the same time units, i.e. Julian days, NSSDC epoch
; serin0  input time series, series of scalars, vectors or  matrices of serin0 = serin0[ .,.,.,nt] nt is the number of timesteps in series
; fillval - fill value
;output
; tout - centered times of binned data ; tout will be in the same time units as tin.
; serout - averaged binned timeseries
; seroutuncertainty - standard deviaton in each bin
; missing_variance_flag - value to indicate a missing sample variance. RY (08/08/2019)
; serout_flag=serout_flag -number of samples per bin

a = serin0
;print,'min/max of a', minmax(a,fill=fillval) ;diagnostic
tin0 = tin00
npts = n_elements(tin0)
;help, fillval, validmin, validmax
cmb_filter_time_series,a,fillval=fillval, validmax=validmax, validmin=validmin ; set values outside of valdimin/max to fillval
;print, 'cmb_timebin_array' & help, a,fillval,validmin,validmax
ii=where( finite(a) eq 0 ) & if ii[0] ne -1 then a[ii]=fillval

if keyword_set(autobad) then begin
   a = cmb_collapse(a, ndimsorg = ndimsorg0)
   ig = cmb_autobad(a,sigmul,fill=fillval, /filter_by_fill)
   ii=where(ig eq 0)
   if ii[0] ne -1 then a[ii]=fillval
   a = cmb_collapse(a, ndimsorg = ndimsorg0,/revert)
endif

if dtbin le 0 then begin ;don't bin data, return full resolution
   tout = tin0
   serout = a
   serout_flag = ''
   return
endif

dt = dtbin
tout = cmb_timesofbin(t0,tf,dt)
ntb = n_elements(tout)

; below if added SAB 6/19/2017
if npts eq 1 then BEGIN ; Ron Yurow pointed out an error is only 1 time instance in timeseries
    print,'***** WARNING ONLY 1 TIME STEP ******'
	serout = serin0 + intarr(n_elements(tout) )
	serout_flag= serout*0
	i = where( tin0 ge tout[0:ntb-2] and tin0 lt tout[1:*])
	if i[0] ne -1 then  serout_flag[i] = 1
	Return 
endif


si0 = size(a,/structure)
if si0.n_dimensions ne 1 then begin ;matrix scalar timeseries
    ndin = si0.dimensions[0:si0.n_dimensions-1]
    nt = ndin[si0.n_dimensions-1] ;time steps
    na = 1l & for i=0,si0.n_dimensions-2 do na = na*ndin[i] ;number of elements in matrix at a timestep
    nd1 = [na,nt]
    a = reform(a,nd1,/overwrite) ; collaspe matrix times series into vector series
    sia = size(a,/structure)
    nc = sia.dimensions[0]
    ndout = ndin & ndout[si0.n_dimensions-1] = ntb ; form internal dimensions of output array
    b = dblarr(na,ntb)*a[0] ;force variable type of b to be the same as a.
endif else begin ;scalar time series
    a = transpose(a)
    b = transpose(dblarr(ntb)*a[0])
    ndout = ntb
    nc=1
endelse
b_nbin = b
b2 = b ;*1d0 ;SAB 9/25/2017 compute uncertainty
iabmap = floor((tin0-t0)/dt)

bin_variance = b

;convert matrix time series into vector time series
for ic = 0l, nc-1 do begin
  ;ig = reform(where( abs(a[ic,*]) lt 0.99*abs(fillval)))
  ig = reform(where( abs(a[ic,*]) ne abs(fillval)))
  ;help,ig & stop
  if ig[0] ne -1 then begin
      ;  This variables are used for the calculation of the variance using 
      ; Welford method
      ; Ron Yurow (July 24, 2018)
      m =    DBLARR (ntb)
      s =    DBLARR (ntb)
      prev = DBLARR (ntb)

      for jpt=0l,n_elements(ig)-1 do begin
		ipt = ig[jpt]
		if iabmap[ipt] ge 0 and iabmap[ipt] lt ntb then begin
		b[ic,iabmap[ipt]]=  b[ic,iabmap[ipt]] + a[ic,ipt]
        ; This method is no longer used for computing variance
        ; Ron Yurow (July 24, 2018)
		; b2[ic,iabmap[ipt]]=  b2[ic,iabmap[ipt]] + a[ic,ipt]^2  ;SAB 9/25/2017 compute uncertainty
		b_nbin[ic,iabmap[ipt]] = b_nbin[ic,iabmap[ipt]] + 1

        ; Use the Welford method to compute the variance.  The Welford method is described in the
        ; following pseudocode
        ; variance(samples):
        ; M := 0
        ; S := 0
        ; for k from 1 to N:
        ;   x := samples[k]
        ;   oldM := M
        ;   M := M + (x-M)/k
        ;   S := S + (x-M)*(x-oldM)
        ; return S/(N-1)
        ; The inner part of the loop is impmlemented below.
        ; Ron Yurow (July 24, 2018)
        prev [iabmap[ipt]] = m [iabmap[ipt]]
        m [iabmap[ipt]] = m [iabmap[ipt]] + (a[ic,ipt] - m [iabmap[ipt]]) / b_nbin [ic,iabmap[ipt]]
        s [iabmap[ipt]] = s [iabmap[ipt]] + (a[ic,ipt] - m [iabmap[ipt]]) * (a [ic,ipt] - prev [iabmap[ipt]])       
		endif
      endfor

      ; Do the actual variance calculation based on the terms previously calculated.
      ; Ron Yurow (July 24, 2018)
      nz = WHERE (b_nbin [ic, *] gt 0)
      bin_variance [ic, nz] = s [nz] / (reform (b_nbin[ic, *]) - 1) [nz]
  endif 
;  if ic eq 5 then stop
endfor

;average
ii=where(b_nbin ne 0)
if ii[0] ne -1 then b[ii] = b[ii]/b_nbin[ii] ; compute mean of each bin


; compute uncertainty of the mean
ii=where(b_nbin gt 2) ;
;if ii[0] ne -1 then b2[ii] = $
;   sqrt((b2[ii]/b_nbin[ii] - b[ii]^2) $
;   *(b_nbin[ii]/(b_nbin[ii]-1))/b_nbin[ii]) ; uncertainty of the mean  ;SAB 9/25/2017 compute uncertainty
if ii[0] ne -1 then BEGIN
    ; No longer needeed, as variance is now calculated using the Welford method.
    ; Ron Yurow (July 24, 2018)
	; sigma2 = (b2[ii] - b_nbin[ii]*b[ii]^2)/(b_nbin[ii]-1) ; standard deviation
	; b2[ii] = sqrt( sigma2 ) ;/b_nbin[ii]) ; comment out uncertainty of mean computatin SAB 3/28/2018
    b2 [ii] = sqrt (bin_variance [ii])  ; Calculate standard deviation based on variance.
endif   
ii=where(b_nbin le 2)
; Changed by Ron Yurow (08/08/2019)
; if ii[0] ne -1 then b2[ii] = fillval ;max(b2) changed  SAB 3/28/2018
if ii[0] ne -1 then b2[ii] = missing_variance_flag
; end compute uncertainty of the mean

;fill in empty bins
if keyword_set(fill_empty_bins) then cmb_fill_array,b,tout,b_nbin else cmb_setemptybin2fill,b,b_nbin,fillval

serout = reform(b,ndout,/overwrite) 
seroutuncertainty = reform(b2,ndout,/overwrite)
 
serout_flag = reform(b_nbin,ndout,/overwrite)

if (where(serout_flag eq 0))(0) eq -1 then emptybin=0 else emptybin = n_elements(where(serout_flag eq 0))
print,'Fraction of missing binned points =  ', float(emptybin)/n_elements(serout_flag)
end