; docformat = 'rst'
;+
; :Author: Paulo Penteado (pp.penteado@gmail.com), Feb/2010
;-
;+
; :Description:
; Examples to test and illustrate the different behaviors of pp_resample.
;
; Makes up a high resolution Gaussian line and resamples at nr locations (defaults to 15).
;
; :Params:
; nr : in, optional, default=15
; The number of locations where the function is to be resampled at.
;
; :Examples:
; Run the pp_resample_test routine, to compare the different methods, at different resampling resolution::
;
; .com pp_resample ;pp_resample_test is in pp_resample.pro, so it will not be found by itself
; pp_resample_test,15 ;At this resolution, the local mode is still reasonable.
; sigma: 7.6632088e-09 6.6129373e-11 1.1358731e-13
; Trapezoidal areas
; New grid Original areas Resampled areas
; -3.700000E+00 0.000000E+00 0.000000E+00
; -3.128571E+00 8.421810E-05 8.421810E-05
; -2.557143E+00 2.646517E-03 2.731632E-03
; -1.985714E+00 4.414658E-02 4.425478E-02
; -1.414286E+00 4.031373E-01 4.035549E-01
; -8.428571E-01 2.067294E+00 2.069880E+00
; -2.714286E-01 6.213192E+00 6.223648E+00
; 3.000000E-01 1.177465E+01 1.179997E+01
; 8.714286E-01 1.579429E+01 1.583423E+01
; 1.442857E+00 1.735853E+01 1.740596E+01
; 2.014286E+00 1.768562E+01 1.773506E+01
; 2.585714E+00 1.772227E+01 1.777198E+01
; 3.157143E+00 1.772447E+01 1.777420E+01
; 3.728571E+00 1.772454E+01 1.777427E+01
; 4.300000E+00 1.772454E+01 1.777427E+01
; Rectangular areas
; Start End Original areas Resampled areas
; -3.985714E+00 -3.414286E+00 1.217446E-05 1.217446E-05
; -3.414286E+00 -2.842857E+00 5.196621E-04 1.176678E-02
; -2.842857E+00 -2.271429E+00 1.176170E-02 1.446414E-01
; -2.271429E+00 -1.700000E+00 1.445446E-01 9.845714E-01
; -1.700000E+00 -1.128571E+00 9.835830E-01 3.834577E+00
; -1.128571E+00 -5.571429E-01 3.829106E+00 9.038319E+00
; -5.571429E-01 1.428571E-02 9.021113E+00 1.415784E+01
; 1.428571E-02 5.857143E-01 1.412464E+01 1.687159E+01
; 5.857143E-01 1.157143E+00 1.682703E+01 1.764558E+01
; 1.157143E+00 1.728571E+00 1.759682E+01 1.776405E+01
; 1.728571E+00 2.300000E+00 1.771449E+01 1.777375E+01
; 2.300000E+00 2.871429E+00 1.772411E+01 1.777418E+01
; 2.871429E+00 3.442857E+00 1.772453E+01 1.777419E+01
; 3.442857E+00 4.014286E+00 1.772454E+01 1.777419E+01
; 4.014286E+00 4.585714E+00 1.772454E+01 1.777419E+01
; pp_resample_test,5 ;At this resolution, the local mode gives an ugly result, but non local is still good.
; sigma: 7.7028521e-09 6.6471473e-11 1.1417492e-13
; Trapezoidal areas
; New grid Original areas Resampled areas
; -3.700000E+00 0.000000E+00 -0.000000E+00
; -1.700000E+00 1.436523E-01 1.436523E-01
; 3.000000E-01 1.177465E+01 1.191919E+01
; 2.300000E+00 1.771441E+01 1.787357E+01
; 4.300000E+00 1.772454E+01 1.788378E+01
; Rectangular areas
; Start End Original areas Resampled areas
; -4.700000E+00 -2.700000E+00 1.201453E-03 1.201453E-03
; -2.700000E+00 -7.000000E-01 2.865225E+00 1.715248E+01
; -7.000000E-01 1.300000E+00 1.714266E+01 1.773729E+01
; 1.300000E+00 3.300000E+00 1.772451E+01 1.773732E+01
; 3.300000E+00 5.300000E+00 1.772454E+01 1.773735E+01
;
; .. image:: pp_resample_test_1.png
; .. image:: pp_resample_test_2.png
;
; :Uses: pp_resample, pp_integral
;
; :Author: Paulo Penteado (pp.penteado@gmail.com), Feb/2010
;-
pro pp_resample_test,nr
compile_opt idl2,hidden
;Defaults
nr=n_elements(nr) eq 1 ? nr : 15
;Make up a well sampled test function and plot it
;This function has so many points that if nr is much smaller than nx, it does not matter whether it is seen
;as a literal or sampled function
nx=10001
x=(dindgen(nx)/(nx-1d0))*16d0-8d0
y=1d1*exp(-x^2)
iplot,x,y,/isotropic,thick=3.,name='Original function (literal)',/insert_legend,$
title='Original function at '+strtrim(string(nx),2)+' points, resampled at '+strtrim(string(nr),2)+' points'
;Make up a grid to resample to
xout=(dindgen(nr)/(nr-1d0))*8d0-3.7d0
areas1=pp_integral(x,y,xmin=replicate(xout[0],nr),xmax=xout,/local)
;Resample with rectangles (local=0) and plot it
yout=pp_resample(y,x,xout,xstart=xstart,xend=xend,xoutbox=xoutbox,youtbox=youtbox)
areas0=pp_integral(x,y,xmin=replicate(xstart[0],nr),xmax=xend)
areas2=pp_integral(xout,yout,/cumulative)
iplot,xoutbox,youtbox,/over,color=[0,0,255],thick=2.,name='Resampled function (sampled)',/insert_legend
iplot,xout,yout,/over,color=[0,255,0],thick=2.,name='Resampled locations (sampled)',/insert_legend,sym_ind=6,sym_size=2.
;Resample with trapezoids (local=1) and plot it
yout2=pp_resample(y,x,xout,/local)
areas3=pp_integral(xout,yout2,/cumulative,/local)
iplot,xout,yout2,/over,color=[255,0,0],name='Resampled function (literal)',/insert_legend,sym_ind=4,sym_size=2.
print,'Trapezoidal areas'
print,'New grid','Original areas','Resampled areas',format='(3A17)'
for i=0,nr-1 do print,xout[i],areas1[i],areas3[i],format='(3E17.6)'
print,'Rectangular areas'
print,'Start','End','Original areas','Resampled areas',format='(4A17)'
for i=0,nr-1 do print,xstart[i],xend[i],areas0[i],areas2[i],format='(4E17.6)'
end
;+
; :Description:
; Resamples the given y(x) function, preserving its area (not just a simple interpolation
; of the function points, as interpol would do). By default (local=0), each input y(x) is
; interpreted a measured sample, that is, as the average of the "flux" falling inside the region
; centered at the corresponding x value: x locations are interpreted as the centers of bins. The same
; interpretation is used for the output locations (xout), and the returned results. See the plots in
; the example for a graphical demonstration.
;
; A literal interpretation of the function (setting the keyword local), as pp_integral does, is intended,
; but not yet properly implemented, so the local mode should not be used at this time.
; The input values must be ordered in x (it does not matter whether it is increasing or decreasing order).
;
; :Returns:
; An array with the same length as xout, with the input function resampled at those locations.
;
; :Params:
; y : in, required
; An array with the function values corresponding to the locations in x.
; x : in, required
; An array of locations where the function is sampled. Must be ordered (increasing or decreasing).
; xout : in, required
; An array of the locations where the y(x) function is to be resampled.
;
; :Keywords:
; xstart : out, optional
; Used when not in local mode, to return the start locations of the bins with center in xout.
; xend : out, optional
; Used when not in local mode, to return the end locations of the bins with center in xout.
; xbox : out, optional
; Used when not in local mode, to return the input x values in pairs of equal values, so that
; plotting xbox,ybox shows the rectangles that are how the input values were interpreted.
; ybox : out, optional
; Used when not in local mode, to return the input y values in pairs of equal values, so that
; plotting xbox,ybox shows the rectangles that are how the input values were interpreted.
; xoutbox : out, optional
; Used when not in local mode, to return the resampled x values in pairs of equal values, so that
; plotting xoutbox,youtbox shows the rectangles that are how the resampled values are to be interpreted.
; youtbox : out,optional
; Used when not in local mode, to return the resampled y values in pairs of equal values, so that
; plotting xoutbox,youtbox shows the rectangles that are how the resampled values are to be interpreted.
; local : in, optional, default=0
; If not set, the y values are interpreted as a measured sample, that is, as the average of the "flux" falling inside the region
; centered at the corresponding x values: x locations are interpreted as the centers of bins. Since this interpretation is equivalent
; to a literal interpretation of a function made of rectangular regions, in this method there are no interpolations (partial areas, if
; any, are just the corresponding fractions of the rectangles). The returned values are to be interpreted in the same way, that is,
; they represent rectangular bins, which have the same area as the rectangular bins in the corresponding region of the input function.
;
; If set, interprets the function literally. The result is a function where the area between each of its
; points (at xout) is the same as the area of y(x) between these locations.
; newton : in, optional, default=0
; Passed on to pp_integral. If local is set and newton is set, the function integrations are done with int_tabulated, which
; uses a 5 point Newton-Cotes formula. If local is not set, this has no effect.
; p0 : in, default=4
; Parameter for the local mode. It should not be necessary to change it from the default, and it should disappear
; in future implementations. It is the number of points to use for the numeric evaluation of a function's
; minimum, but the result is not expected to vary for p0 larger than 3, as the function is expected to be
; an exact parabola.
; lsquadratic : in, optional, default=0
; Passed on to pp_integral. If the method requires interpolation, this is passed on to interpol,
; to select 4 point quadratic interpolation.
; If none of lsquadratic, quadratic and spline are set, interpolation is linear.
; quadratic : in, optional, default=0
; Passed on to pp_integral. If the method requires interpolation, this is passed on to interpol,
; to select 3 point quadratic interpolation.
; If none of lsquadratic, quadratic and spline are set, interpolation is linear.
; spline : in, optional, default=0
; Passed on to pp_integral. If the method requires interpolation, this is passed on to interpol,
; to select spline interpolation.
; If none of lsquadratic, quadratic and spline are set, interpolation is linear.
;
; :Examples:
;
; Make a well-sampled input function y(x) and an output grid xout with few locations::
;
; x=dindgen(10001)/1d4-0.5d0
; y=exp(-(x*2d1)^2)*0.5d0
; xout=dindgen(5)/6d0-0.3d0
;
; Do the resampling and look at the results::
;
; yout=pp_resample(y,x,xout,xstart=xstart,xend=xend,xbox=xbox,ybox=ybox,xoutbox=xoutbox,youtbox=youtbox)
; print,xstart,xend
; ;-0.38333333 -0.21666667 -0.050000000 0.11666667 0.28333333
; ;-0.21666667 -0.050000000 0.11666667 0.28333333 0.45000000
; print,xstart,xend
; ;-0.38333333 -0.21666667 -0.050000000 0.11666667 0.28333333
; ;-0.21666667 -0.050000000 0.11666667 0.28333333 0.45000000
; iplot,xbox,ybox,name='input function (rectangles)',/insert_legend,thick=2.
; iplot,xoutbox,youtbox,/over,color=[255,0,0],name='resampled function (rectangles)',/insert_legend,thick=2.
; iplot,xout,yout,/over,color=[0,255,0],name='resampled function (locations)',/insert_legend,thick=2.
;
; .. image:: pp_resample.png
;
; See also the example in pp_resample_test, for a comparison of the methods.
;
; :Uses: pp_integral
;
; :Todo:
; Improve the implementation of the resampling when local is set, to get replace the numerical
; minimization by an analytic expression.
;
;
; :Author: Paulo Penteado (pp.penteado@gmail.com), Feb/2010
;-
function pp_resample,y,x,xout,xstart=xstart,xend=xend,xbox=xbox,ybox=ybox,xoutbox=xoutbox,youtbox=youtbox,$
local=local,newton=newton,p0=p0,$
lsquadratic=lsquadratic,quadratic=quadratic,spline=spline
compile_opt idl2
;Check that the input is valid
nx=n_elements(xout)
if (nx lt 2L) then message,'Output grid xout must contain more than two points'
;Defaults
local=n_elements(local) eq 1 ? local : 0
p0=n_elements(p0) eq 1 ? p0 : 4
if (~local) then begin ;Interpret the functions as sample averages
;Calculate the start and en of each output point
xstart=[xout[0]-(xout[1]-xout[0])/2d0,(xout[0:nx-2]+xout[1:nx-1])/2d0]
xend=[xstart[1:nx-1],xout[nx-1]+(xout[nx-1]-xout[nx-2])/2d0]
;Get area between the start and end of each output point
areas=pp_integral(x,y,xmin=xstart,xmax=xend,xstart=ixstart,xend=ixend,yinc=yinc,$
local=local,newton=newton,lsquadratic=lsquadratc,quadratic=quadratic,spline=spline)
;Calculate the function average over the span of each output point
ret=areas/(xend-xstart)
;Create xoutbox and youtbox arrays, with constant y values spanning the start and end of each xout point
if (arg_present(xoutbox) || arg_present(youtbox)) then begin
xoutbox=dblarr(2L*nx) & youtbox=dblarr(2L*nx)
xoutbox[0:*:2]=xstart & xoutbox[1:*:2]=xend
youtbox[0:*:2]=ret & youtbox[1:*:2]=ret
endif
endif else begin ;Interpret the functions literally
;Calculate the cumulative areas at each point on the input function
areas=pp_integral(x,y,/local,xmin=xout[0:nx-2],xmax=xout[1:nx-1],newton=newton,lsquadratic=lsquadratc,quadratic=quadratic,spline=spline)
;Numerical determination of the smallest line with the same area as the input function
yf=sqrt(total(y^2))
yf=yf gt 0d0 ? yf : 1d0
r0=yf*dindgen(p0)/(p0-1d0)*2d0-1d0
tmp=2d0*areas[0:nx-2]/(xout[1:nx-1]-xout[0:nx-2])
ret=dblarr(nx)
rms=dblarr(p0)
for i=0,p0-1 do begin
ret[0]=r0[i]
for j=1,nx-1 do ret[j]=tmp[j-1]-ret[j-1]
rms[i]=total(ret^2)
endfor
;rms(r0) should be a parabola, so fit a parabola to find its minimum
pol=poly_fit(r0,rms,2,/double,sigma=sigma)
print,'sigma: ',sigma ;Uncertainties on the parabola parameters
;Recalculate the line with the r that produces the smallest line
ret[0]=-pol[1]/(2d0*pol[2])
for j=1,nx-1 do ret[j]=tmp[j-1]-ret[j-1]
endelse
;Create xbox and ybox arrays, with constant y values spanning the start and end of each x point
if (arg_present(xbox) || arg_present(ybox)) then begin
nix=n_elements(x)
xbox=dblarr(2L*nix) & ybox=dblarr(2L*nix)
xbox[0:*:2]=ixstart & xbox[1:*:2]=ixend
ybox[0:*:2]=yinc & ybox[1:*:2]=yinc
endif
return,ret
end