function alaloess,yy,xx,hwidth=hwidth,ysig=ysig,ndeg=ndeg,\$ yerr=yerr,verbose=verbose, _extra=e ;+ ;function alaloess ; computes a loess/lowess type smoothed function for the ; given set of input points by fitting polynomials over ; the neighborhood of each point and using the local value ; of the polynomial fit as the smoothed curve ; ; another one in a hoary line of PINTofALE smoothing routines, ; joining SMOOTHIE, NOISMOOTH, VARSMOOTH, VOORSMOOTH, REGROUP, ; SPLAC, CONV_RMF, SCRMF, UNKINK, LINEREM, HAARLINE, etc. ; ;syntax ; yl=alaloess(yy,xx,hwidth=hwidth,ysig=ysig,ndeg=ndeg,\$ ; yerr=yerr,verbose=vv) ; ;parameters ; yy [INPUT; required] the points that need to be smoothed ; xx [INPUT] the abscissae for which YY are defined ; * if size eq size of YY or +1, taken to be the midpoints ; or the grid ; * if insensible, taken to be array indices, LINDGEN(N(YY)) ; ;keywords ; hwidth [INPUT] half range over which local polynomial fits are ; to be performed ; * if XX are array indices, HWIDTH is forced to be >0 ; * if XX is a predefined grid, the default value is the ; minimum of delta(XX) ; ysig [INPUT] error on YY ; * if size matches that of YY, then maps 1-to-1 ; * if scalar, taken to be ; <0 : absolute error ; >0 and <1 : fractional error ; >1 : percentage error ; =0 : ignored in call to POLY_FIT() ; on YY ; * otherwise, assumed to be sqrt(abs(YY)+0.75)+1 ; ndeg [INPUT] degree of the polynomial to be fit ; * default is 3 ; yerr [OUTPUT] 1-sigma error at the smooth location from ; each fit. ; * WARNING: this is to be taken only as an illustrative ; number and not used for anything quantitative, because ; the errors in adjacent bins are not independent ; verbose [INPUT] controls chatter ; _extra [JUNK] here only to prevent crashing the program ; ;example ; .run alaloess ;(requires PEASECOLR) ; ;these are all the smoothing tools in PINTofALE ; ALALOESS() makes a loess curve ; CLTSMOOTH() removes -ves from background subtracted spectra or light curves ; CONV_RMF convolves with specified RMF ; HAARTRAN() smoothes by filtering on Haar wavelet coefficients ; LINEREM() removes lines from spectrum by iterative S/N filtering ; NOISMOOTH() does boxcar accumulation a la Ebeling's asmooth ; REGROUP() accumulates from one end ; SCRMF does fast convolution using optimized RMF ; SMOOTHIE does peak-descent accumulation to build up S/N ; SPLAC computes a segmented piecewise linear approximations ; UNKINK() removes sharp kinks from a curve ; VARSMOOTH() does local smoothing at specified scales ; VOORSMOOTH() does fast local smoothing at specified scales ; ;history ; vinay kashyap (MarMMX) ; some bug fixes, and hardening to poly_fit failures -- where it fails ; now simply keep old value (VK; MayMMX) ; updated to be more robust to NaNs (VK; AprMMV) ;- ; usage np=n_params() & ok='ok' & ny=n_elements(yy) & nx=n_elements(xx) if np eq 0 then ok='Insufficient parameters' else \$ if ny eq 0 then ok='Y(X) is undefined' else \$ if ny eq 1 then ok='Y(X) must be a vector' if ok ne 'ok' then begin print,'Usage: yl=alaloess(y,x,hwidth=hwidth,ysig=ysig,ndeg=ndeg,\$' print,' yerr=yerr,verbose=vv)' print,' returns lowess like smooth curve based on local polynomial fit' if np ne 0 then message,ok,/informational return,-1L endif ; parameters yz=0.0+[yy[*]] & xok='units' case nx of ny: xz=xx ny+1: xz=0.5*(xx[1:*]+xx) else: begin xz=lindgen(ny) xok='indices' endelse endcase ; keywords vv=0L & if keyword_set(verbose) then vv=long(verbose[0])>1L ; mpoly=3 & if keyword_set(ndeg) then mpoly=abs(fix(ndeg[0])) ; nye=n_elements(ysig) & sigy=sqrt(abs(yz)+0.75)+1. if nye eq 1 then begin if ysig[0] lt 0 then sigy[*]=abs(ysig[0]) if ysig[0] gt 0 and ysig[0] lt 1 then sigy=ysig[0]*yz if ysig[0] ge 1 then sigy=(ysig[0]/100.)*yz if ysig[0] eq 0 then sigy[*]=0 endif if nye eq ny then sigy=ysig+0.0 ; hrange=min(abs(xz[1:*]-xz)) if keyword_set(hwidth) then hrange=hrange>hwidth[0] ; step through the array and get the smoothed values yl=yz & yerr=yl if 2*hrange le mpoly then begin message,'fitting window too small for polynomial',/informational message,'returning without any smoothing',/informational print,mpoly,hrange yerr=sigy return,yl endif for i=0L,ny-1L do begin i0=(i-hrange)>0 i1=(i+hrange)<(ny-1L) if xok eq 'units' then begin oo=where(abs(xz-xz[i]) le hrange,moo) i0=oo[0] & i1=oo[moo-1L] endif j=i-i0 if i1-i0 gt mpoly then begin xin=xz[i0:i1] & yin=yz[i0:i1] & yein=sigy[i0:i1] if sigy[i] gt 0 then begin tmp=poly_fit(xin,yin,mpoly,measure_errors=yein,yfit=yfit,\$ yband=yband,/double,status=ist) if ist gt 0 then tmp=poly_fit(xin,yin,mpoly,yfit=yfit,yband=yband,/double,status=ist) endif else begin tmp=poly_fit(xin,yin,mpoly,yfit=yfit,yband=yband,/double,status=ist) if ist gt 0 then tmp=poly_fit(xin,yin,mpoly,/double,status=ist) endelse if ist eq 0 then begin yl[i]=yfit[j] & yerr[i]=yband[j] endif endif endfor if vv gt 1000 then stop,'HALTing; type .CON to continue' return,yl end ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; example calling sequence yl=alaloess() if not keyword_set(ndeg) then ndeg=3 if not keyword_set(hwidth) then hwidth=8 if n_elements(verbose) eq 0 then verbose=1 print,'HWIDTH='+strtrim(hwidth,2) print,'NDEG='+strtrim(ndeg,2) print,'VERBOSE='+strtrim(VERBOSE,2) y=50*abs(randomn(seed))*mk_gauss(findgen(256),128,10,1)+\$ 30*abs(randomn(seed))*mk_gauss(findgen(256),165,15,1) yoff=40*randomu(seed)+10. yr=long(0*y) for i=0,256-1 do yr[i]=randomu(seed,poisson=y[i]+yoff) yl=alaloess(yr,hwidth=hwidth,ndeg=ndeg,yerr=yerr,verbose=verbose) peasecolr & loadct,3 & peasecolr plot,yr,psym=1 & oplot,y+yoff,col=1 oplot,yl,col=2 & oplot,yl+yerr,col=2,line=2 & oplot,yl-yerr,col=2,line=2 print,'' & print,'blue: actual; red: smoothed' & print,'' end