;+
;function	rundens
;	compute and return the running density of a set of measurements.
;
;syntax
;	yrun=rundens(xpt,hwidth,/usex,/optwdt,outhwdt=outhwdt,$
;	yerr=yerr,xerr=xerr,nsim=nsim,ysim=ysim,yesim=yesim,xtol=xtol,$
;	wts=wts,dwts=dwts,lbound=lbound,rbound=rbound,cvmse=cvmse,$
;	verbose=verbose)
;
;parameters
;	xpt	[INPUT; required] array of points for which the running density
;		must be calculated
;		* a running density is calculated at each value of XPT
;	hwidth	[INPUT] the half-width on either side of a point to include in
;		the density
;		* this is assumed to be the number of array indices
;		  unless USEX is set
;		* default is to use 5% of n_elements(XPT)
;		  or 5% of the range of XPT if USEX is set
;		* if OPTWDT is set, computes the optimum width to use
;		  and treats this as just a default
;		* actual value used is returned via keyword OUTHWDT
;
;keywords
;	usex	[INPUT] if set, assumes HWIDTH is in same units as XPT
;	optwdt	[INPUT] if set, figures out the optimum HWIDTH to use by
;		comparing the width of the distribution of the running
;		density estimates to the expected statistical stddev for
;		a flat density
;		* if this is set, HWIDTH will be overridden
;		* the optimum HWIDTH is calculated based on the array
;		  indices -- i.e., USEX is ignored during the calculation,
;		  but if set, is used afterwards to scale the optimum
;		  HWIDTH to the units in XPT.  So beware that it may not
;		  really be the optimum in this case.
;		* if this module fails, the smallest width checked will
;		  be returned
;	outhwdt	[OUTPUT] will contain the value of HWIDTH actually used in
;		the calculation
;	yerr	[OUTPUT] error bar on the density at each point
;		* WARNING: adjacent values of YERR are not independent.
;		  while the error estimate at each point is fine, to
;		  get a truly independent estimate, one has to skip by
;		  2*HWIDTH
;	xerr	[INPUT] error bars on XPT
;		* used to compute a Monte Carlo uncertainty band on the
;		  output running density
;	nsim	[INPUT] number of simulations to carry out
;		* default is none
;	ysim	[OUTPUT] a NSIMxN(XPT) array to contain the running
;		density for each simulation
;	yesim	[OUTPUT] point-wise stddev based on YSIM
;	xtol	[INPUT] the tolerance for considering two points to be
;		identical
;		* essentially, if any two adjacent points are found to
;		  be identical, then all points are jittered by +-0.5*XTOL
;	wts	[INPUT] if given, computes the total of the WTS and returns
;		that rather than the running density
;		* MUST match the size of XPT, otherwise ignored
;	dwts	[INPUT] same as WTS, but the total is computed and placed
;		in the denominator 
;		* MUST match the size of XPT, otherwise ignored
;		* if both WTS and DWTS are given, the ratio of the totals
;		  is computed, not the average of the ratios
;	lbound	[OUTPUT] array indices of leftside bounds for all XPT
;	rbound	[OUTPUT] array indices of rightside bounds for all XPT
;	cvmse	[OUTPUT] cross-validation mean square error between estimate
;		and leave-one-out estimate.  the extra calculations are
;		performed only if this keyword is present and are not done
;		for the simulated data.
;	verbose	[INPUT] controls chatter
;	
;	_extra	[JUNK] here only to prevent crashing the program
;
;subroutines
;	RUNDENS_OPTWDT (included here)
;	CURVESECT
;	EQT_INTERVAL
;
;history
;	vinay kashyap (MMXI.VII)
;	added keywords WTS and DWTS and corrected bug (VK; MMXI.VIII)
;	added keywords LBOUND, RBOUND, and CVMSE (VK; MMXII.V)
;	added keyword OPTWDT, added function RUNDENS_OPTWDT, and added
;	  calls to RUNDENS_OPTWDT(), CURVESECT(), and EQT_INTERVAL()
;	  (VK; MMXII.XII)
;-

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
function rundens_optwdt,xpt,verbose=verbose, _extra=e
;{
;function	rundens_optwdt
;	repeatedly calls rundens() with different values of HWIDTH
;	to find the optimal HWIDTH and returns that value
;
;	the optimum is defined as that value of HWIDTH where the
;	scatter in the running means is comparable to the expected
;	scatter for the given number of points in a flat density curve.
;
;dedicated function, expected to called only from RUNDENS()
;
;vinay kashyap (MMXII.XII)
;}

;	keywords
vv=0L & if keyword_set(verbose) then vv=long(verbose[0])>1L

nph=n_elements(xpt)
if nph eq 0 then begin
  message,'XPT is not defined; quitting',/informational & return,-1L
endif
if nph lt 3 then begin
  message,'XPT is too small; quitting',/informational & return,-1L
endif
if nph lt 16 then begin
  message,'XPT too small; half the distance to the goal line',/informational
  return,nph/2
endif

kph2=fix(alog(nph)/alog(2))
kk=round(sqrt(2)^(findgen(kph2*2-6)+6))
nkk=n_elements(kk)
dfvar=fltarr(nkk) & sfvar=fltarr(nkk)
for k=0L,nkk-1L do begin
  kph=kk[k]
  f=rundens(xpt,kph, _extra=e)
  df=f[1:*]-f
  ;dfvar[k]=stddev(df)	;this is too conservative
  tmp=eqt_interval(df,/fsample,clev=0.68) & dfvar[k]=(tmp[1]-tmp[0])/2.
  sfvar[k]=stddev(f)/sqrt(nph)
endfor

xmax=max(xpt,min=xmin) & dx=xmax-xmin
tmp=sqrt(2*kk+1L)/dx & tmp=sqrt(tmp^2+sfvar^2)
;xy=curvesect(dfvar,sqrt(2*kk+1L)/dx,kk)
xy=curvesect(dfvar,tmp,kk)
if finite(xy[1]) eq 0 then begin
  if vv gt 0 then message,'No intersections found; trying conservative match',$
  	/informational
  tmp=sqrt(2*kk+1L)/dx
  xy=curvesect(dfvar,tmp,kk)
  if finite(xy[1]) eq 0 then begin
    message,'No intersections found; choosing the smallest value',$
    	/informational
    xy[0]=min(kk)
  endif
endif
hw=long(xy[0])

if vv gt 2000 then stop,'HALTing; type .CON to continue'

return,hw
end

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
function rundens,xpt,hwidth,usex=usex,optwdt=optwdt,outhwdt=outhwdt,$
	yerr=yerr,xerr=xerr,nsim=nsim,ysim=ysim,yesim=yesim,xtol=xtol,$
	wts=wts,dwts=dwts,lbound=lbound,rbound=rbound,cvmse=cvmse,$
	verbose=verbose, _extra=e

;	usage
ok='ok' & np=n_params() & nx=n_elements(xpt) & nw=n_elements(hwidth)
if np lt 2 then ok='Insufficient parameters' else $
 if nx eq 0 then ok='XPT is undefined' else $
  if nx eq 1 then ok='XPT is not an array'
if ok ne 'ok' then begin
  print,'Usage: yrun=rundens(xpt,hwidth,/usex,/optwdt,outhwdt=outhwdt,$
  print,'       yerr=yerr,xerr=xerr,nsim=nsim,ysim=ysim,yesim=yesim,xtol=xtol,$
  print,'       wts=wts,dwts=dwts,lbound=lbound,rbound=rbound,cvmse=cvmse,$
  print,'       verbose=verbose)
  print,'  compute and return running density'
  if np ne 0 then message,ok,/informational
  return,-1L
endif

;	inputs
vv=0L & if keyword_set(verbose) then vv=long(verbose[0])>1
;
xmin=min(xpt,max=xmax)
hw=long(0.05*nx+0.5) & if keyword_set(usex) then hw=0.05*(xmax-xmin)
if nw ne 0 then hw=1.*abs(hwidth[0]) & if not keyword_set(usex) then hw=long(hw)
nw=n_elements(hw)
;
msim=0L & if keyword_set(nsim) then msim=long(abs(nsim[0]))>1
;
xsig=fltarr(nx) & nxe=n_elements(xerr)
if nxe gt 0 then xsig[0L:nxe<(nx-1L)]=abs(xerr[0L:nxe<(nx-1L)])
;
tolx=1d-10 & if keyword_set(xtol) then begin
  if xtol[0] ne 1 then tolx=abs(xtol[0])
endif
;
numwt=fltarr(nw)+1. & denwt=numwt
nnwt=n_elements(wts) & if nnwt ne nx then nnwt=0L else numwt=1.0*wts
ndwt=n_elements(dwts) & if ndwt ne nx then ndwt=0L else denwt=1.0*dwts
;
if keyword_set(optwdt) then begin
  ;opthw=rundens_optwdt(xpt,xtol=tolx,wts=numwt,dwts=denwt,verbose=vv) ;UNCOMMENT THIS ONLY AFTER IMPLEMENTING YERR FUNCTIONALITY IN RUNDENS_OPT
  opthw=rundens_optwdt(xpt,xtol=tolx,verbose=vv)
  if opthw[0] ne -1 then begin
    if keyword_set(usex) then hw=opthw*(xmax-xmin)/float(nx) else hw=opthw
    if vv gt 5 then message,$
    	'using HWIDTH='+strtrim(hwidth,2),/informational
  endif else if vv gt 0 then message,$
    	'OPTWDT did not work; continuing with inputs and/or defaults',$
	/informational
endif

;	outputs
outhwdt=hw
yrun=fltarr(nx)*xpt[0] & yerr=yrun
if msim gt 0 then begin
  ysim=xpt[0]*fltarr(msim,nx) & yesim=xpt[0]*fltarr(nx)
endif
lbound=lonarr(nx)-1L & rbound=lbound
if arg_present(cvmse) then begin
  crossval=1
  yrun2=fltarr(nx)*xpt[0] & yerr2=yrun2
endif else crossval=0

;	step through the array and compute the running density
for isim=0L,msim do begin	;{isim=0 is standard run, rest are sims

  xx=xpt & if isim gt 0 then xx=xx+randomn(seed,nx)*xsig
  zn=numwt & zd=denwt
  os=sort(xx) & xx=xx[os] & zn=zn[os] & zd=zd[os]
  delx=min(xx[1:*]-xx)
  if delx eq 0 then begin
    xx=xx+randomu(seed)*(tolx/2.)-tolx/2.
    os=sort(xx) & xx=xx[os] & zn=zn[os] & zd=zd[os]
  endif
  if isim eq ((100-vv)>1)*long(isim/((100-vv)>1)) then kilroy

  for i=0L,nx-1L do begin		;{step through the input array

    if keyword_set(usex) then begin	;(find array range over which to compute density
      ok=where(abs(xx-xx[i]) le hw,mok)
      i0=ok[0] & i1=ok[mok-1L]
      if mok gt 0 then begin
        yy=float(mok)
        if isim eq 0 then ye=sqrt(mok) & yne=0 & yde=0
        ;
        if nnwt gt 0 then begin & yn=total(zn[ok]) & if mok gt 1 and isim eq 0 then yne=stddev(zn[ok]) & endif
        if ndwt gt 0 then begin & yd=total(zd[ok]) & if mok gt 1 and isim eq 0 then yde=stddev(zd[ok]) & endif
        if nnwt gt 0 and ndwt eq 0 then begin
          yy=yn & if isim eq 0 then ye=sqrt(yne^2+mok*(yy/mok)^2)	;E(var)+var(E), with E(var)~stddev and var(E)~N*<y>^2
        endif
        if nnwt eq 0 and ndwt gt 0 then begin
          yy=1./yd & if isim eq 0 then ye=yy*sqrt(yde^2/yd^4+mok/(yy/mok)^2)
        endif
        if nnwt gt 0 and ndwt gt 0 then begin
          yy=yn/yd & if isim eq 0 then ye=yy*sqrt((yne/yd)^2+(yn*yde/yd^2)^2+mok*(yy/mok)^2)
        endif
        ;
        yy=yy/hw/2.
        if isim eq 0 then ye=ye/hw/2.
	;
	if isim eq 0 and keyword_set(crossval) then begin	;(compute array for cross-validation
	  ok0=where((xx-xx[i]) ge -hw and (xx-xx[i]) le 0,mok0)
	  ok1=where((xx-xx[i]) ge 2*hw and (xx-xx[i]) le 3*hw,mok1)
	  j00=ok0[0] & j01=ok0[mok0-1L] & j10=ok1[0] & j1=ok1[mok1-1L]
	  if mok0 gt 0 and mok1 gt 0 then ok2=[ok0,ok1] else $
	   if mok0 gt 0 and mok1 eq 0 then ok2=ok0 else $
	    if mok0 eq 0 and mok1 gt 0 then ok2=ok1 else $
	     ok2=[-1L] & mok2=mok0+mok1
	  if mok2 gt 0 then begin
	    yy2=float(mok2)
	    ye2=sqrt(mok2) & yne2=0 & yde2=0
	    ;
	    if nnwt gt 0 then begin & yn2=total(zn[ok2]) & if mok2 gt 1 then yne2=stddev(zn[ok2]) & endif
	    if ndwt gt 0 then begin & yd2=total(zd[ok2]) & if mok2 gt 1 then yde2=stddev(zd[ok2]) & endif
	    if nnwt gt 0 and ndwt eq 0 then begin
	      yy2=yn2 & ye2=sqrt(yne2^2*mok2*(yy2/mok2)^2)
	    endif
            if nnwt eq 0 and ndwt gt 0 then begin
              yy2=1./yd2 & ye2=sqrt(yde2^2/yd2^4+mok2/(yy2/mok2)^2)
            endif
            if nnwt gt 0 and ndwt gt 0 then begin
              yy2=yn2/yd2 & ye2=sqrt((yne2/yd2)^2+(yn2*yde2/yd2^2)^2+mok2*(yy2/mok2)^2)
            endif
	    ;
	    yy2=yy2/hw/2.
	    ye2=ye2/hw/2.
	  endif
	endif							;ISIM=0 && CROSSVAL1)
      endif
    endif else begin			;USEX)(do not USE X
      i0=(i-hw)>0
      i1=(i+hw)<(nx-1L)
      mok=i1-i0+1L
      dx=xx[i1]-xx[i0]
      yy=float(mok) & if isim eq 0 then ye=sqrt(mok) & yne=0 & yde=0
      ;
      if nnwt gt 0 then begin & yn=total(zn[i0:i1]) & if mok gt 1 and isim eq 0 then yne=stddev(zn[i0:i1]) & endif
      if ndwt gt 0 then begin & yd=total(zd[i0:i1]) & if mok gt 1 and isim eq 0 then yde=stddev(zd[i0:i1]) & endif
      if nnwt gt 0 and ndwt eq 0 then begin
        yy=yn & if isim eq 0 then ye=sqrt(yne^2+mok*(yy/mok)^2)
      endif
      if nnwt eq 0 and ndwt gt 0 then begin
        yy=1./yd & if isim eq 0 then ye=yy*sqrt(yde^2/yd^4+mok/(yy/mok)^2)
      endif
      if nnwt gt 0 and ndwt gt 0 then begin
        yy=yn/yd & if isim eq 0 then ye=yy*sqrt((yne/yd)^2+(yn*yde/yd^2)^2+mok*(yy/mok)^2)
      endif
      ;
      if dx gt 0 then yy=yy/dx
      if dx gt 0 and isim eq 0 then ye=ye/dx
      ;
      if isim eq 0 and keyword_set(crossval) then begin	;(compute array for cross-validation
        j00=i0 & j01=i & j10=(i+2*hw)<(nx-1) & j11=(i+3*hw)<(nx-1)
	mok0=j01-j00+1L & mok1=j11-j10+1L & mok2=mok0+mok1
	dx0=xx[j01]-xx[j00] & dx1=xx[j11]-xx[j10] & dx2=dx0+dx1
	yy2=float(mok2) & ye2=sqrt(mok2) & yne2=0 & yde2=0
	;
	if mok0 gt 0 and mok1 gt 0 then zn2=[zn[j00:j01],zn[j10:j11]] else $
	 if mok0 gt 0 and mok1 eq 0 then zn2=zn[j00:j01] else $
	  if mok0 eq 0 and mok1 gt 0 then zn2=zn[j10:j11]
	if mok0 gt 0 and mok1 gt 0 then zd2=[zd[j00:j01],zd[j10:j11]] else $
	 if mok0 gt 0 and mok1 eq 0 then zd2=zd[j00:j01] else $
	  if mok0 eq 0 and mok1 gt 0 then zd2=zd[j10:j11]
        if nnwt gt 0 then begin & yn2=total(zn2) & if mok2 gt 1 then yne2=stddev(zn2) & endif
        if ndwt gt 0 then begin & yd2=total(zd2) & if mok2 gt 1 then yde2=stddev(zd2) & endif
        if nnwt gt 0 and ndwt eq 0 then begin
          yy2=yn2 & ye2=sqrt(yne2^2+mok2*(yy2/mok2)^2)
        endif
        if nnwt eq 0 and ndwt gt 0 then begin
          yy2=1./yd2 & ye2=yy2*sqrt(yde2^2/yd2^4+mok2/(yy2/mok2)^2)
        endif
        if nnwt gt 0 and ndwt gt 0 then begin
          yy2=yn2/yd2 & ye2=yy2*sqrt((yne2/yn2)^2+(yn2*yde2/yd2^2)^2+mok2*(yy2/mok2)^2)
        endif
        ;
        if dx gt 0 then yy2=yy2/dx2
        if dx gt 0 then ye2=ye2/dx2
      endif						;ISIM=0 && CROSSVAL)
    endelse				;I0,I1)
    lbound[i]=i0 & rbound[i]=i1

    if isim eq 0 then begin
      yrun[i]=yy
      yerr[i]=ye
      if keyword_set(crossval) then yrun2[i]=yy2
      if keyword_set(crossval) then yerr2[i]=ye2
    endif else ysim[isim-1L,i]=yy

  endfor				;I=0,NX-1}
endfor				;ISIM=0,MSIM}
if msim gt 0 then for i=0L,nx-1L do yesim[i]=stddev(ysim[*,i])

;	now compute the cross-validation mean-square error
;	for the running mean, there are K=NX/(2*HW) independent points.
;	so pick these K points and compute mean of squared difference
;	between estimate vs gapped estimate.  The gapped estimates are
;	simply the ones where the Kth bin is removed and the running
;	mean is computed by conflating the next bin with the previous one.
if keyword_set(crossval) then begin
  if keyword_set(usex) then kk=long((max(xx)-min(xx))/hw/2.+1)>1 else kk=(1+nx/hw/2)>1
  jkk=lindgen(kk+1)*float(nx)/float(kk+1)
  ;ikk=[(jkk-1L)>0,jkk,(jkk+1L)<(nx-1L)] & ikk=ikk[uniq(ikk,sort(ikk))]
  ikk=jkk
  yk=yrun[ikk] & yk2=yrun2[ikk] & yek=yerr[ikk]
  ;cvmse=mean(((yk-yk2)/yek)^2,/nan)	;don't do this, it'll "correct" for a crucial factor!
  cvmse=mean(((yk-yk2))^2,/nan)
endif

if vv gt 1000 then stop,'HALTing; type .CON to continue'

return,yrun
end