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