FUNCTION interp_pruned_x, x, ys, ac, ABS_ERROR = abs_error, $
PLOT_PRUNE=plot_prune
;
; This procedure is used to "prune" a (possibly multi-Y-column)
; lookup table by returning the indices of the xs array which,
; if they alone are extracted from the table and used with
; linear interpolation, will lead to interpolated values with
; errors less than some specified level (for each Y-column).
;
; EXAMPLE of USE
;
; hetgcal> meg = rdb_read(!DDHETGCAL+'/cip/hetgmegD1996-11-01greffN0002.rdb')
; hetgcal> acs = [0.001, 0.003, 0.01, 0.03, 0.1]
; hetgcal> for ia=0,4 do xs_needed = interp_pruned_x(meg.energy, $
; [[meg.op1], [meg.oz]], acs(ia)*[1.,1.])
;
; Need 791 points out of 5001 to get accuracy of 0.100000 %
;
; Need 441 points out of 5001 to get accuracy of 0.300000 %
;
; Need 250 points out of 5001 to get accuracy of 1.00000 %
;
; Need 142 points out of 5001 to get accuracy of 3.00000 %
;
; Need 74 points out of 5001 to get accuracy of 10.0000 %
;
;
; INPUTS
;
; x - the "X" array of the lookup table
;
; ys - One or more arrays that correspond to looked-up values
; w.r.t. x, that is ys(*,0), ..., y(*,N-1) are each arrays
; with the same lenght as x. ys is conveniently formed as:
; ys = [ [lookup1], [lookup2], ..., [lookupN-1] ]
;
; as - are the desired accuracy for each of the y's in relative
; terms, that is:
;
; | (pruned_Y(m) - actual_Y(m) ) / actual_Y(m) | < as(m)
;
; RETURNED VALUE
; The indices of the rows in the table that are needed to achieve the
; desired accuracy are returned.
; KEYWORDS
;
; ABS_ERROR - if this keyword is set then the as values are
; interpreted as error values in an absolute sense,
; that is:
;
; | (pruned_Y(m) - actual_Y(m) ) | < as(m)
;
; FUNCTION REQUIRED
; interp_prune_err(ixl, ixn, x, ys) - calculates and returns the
; max error for each table when interpolated from indices ixl to ixn.
;
; ALGORITHM overview: ixl points at the index of that last point
; which we've included in the pruned x list and ixn points at
; the candidate next point. The errors for the tables if
; linear interpolation is used between ixl and ixn is calculated.
; If all of the tables have errors less than desired then
; the value of ixn is incremented by one and we repeat, otherwise
; (errors are too big in at least one table) we add the point
; ixn-1 to the pruned-x-list and set ixl = ixn-1 and repeat.
;
; --------------------------------------------------------------
; 1999-01-16 dd Initial version created.
; --------------------------------------------------------------
;
; Make sure these values are set so they can be passed...
if NOT(KEYWORD_SET(ABS_ERROR)) then abs_error = 0
if NOT(KEYWORD_SET(PLOT_PRUNE)) then plot_prune = 0
nxs = n_elements(x)
; Start by setting ixl to 0 and xs_needed to [0]:
xs_needed = [0]
ixl=0
; Here is the general loop to keep increasing ixn until the
; error is too large, then set ixl to that ixn-1 and continue...
; until ixl is equal to the last x value (always accepted).
while ixl LT nxs-1 do begin
; Show the progress...
if KEYWORD_SET(PLOT_PRUNE) then print, $
' '+STRCOMPRESS(ixl)+'/'+STRCOMPRESS(nxs)
; catch the case that ixl is the next-to-last point, if so
; accept it... and next iteration we'll be done
if ixl EQ nxs-2 then begin
; add the last point
xs_needed = [xs_needed, (ixl+1)]
; and move ixl:
ixl = ixl+1
end else begin
; OK, now we can set ixn to ixl+1
ixn = ixl + 1
; and start increasing ixn as long as the errors are OK
err_OK = 1
while err_OK do begin
ixn = ixn + 1
; if ixn points outside the array then we're done -
; signal an error and ixl will be set to ixn-1
if ixn GE nxs then begin
err_OK = 0
end else begin
max_errors = interp_prune_err(ixl, ixn, x, ys, $
PLOT_PRUNE=plot_prune, ABS_ERROR=abs_error)
err_OK = ( MAX(max_errors GT ac) EQ 0 )
end
end
; at this point ixn points at a point which gives too
; large an error, so set ixl = ixn-1 and continue...
ixl = ixn-1
xs_needed = [xs_needed, (ixl)]
end
end ; of going through the x's.
if KEYWORD_SET(ABS_ERROR) then begin
print, ' Need ' + STRCOMPRESS(n_elements(xs_needed)) + ' points out of ' + $
STRCOMPRESS(nxs)+ ' to get accuracy of '+STRCOMPRESS(ac(0))
end else begin
print, ' Need ' + STRCOMPRESS(n_elements(xs_needed)) + ' points out of ' + $
STRCOMPRESS(nxs)+ ' to get accuracy of '+STRCOMPRESS(100.*ac(0))+' %'
end
print, ''
RETURN, xs_needed
END