spd_ui_interpolate.pro (Documentation for /home/davin/idl/socware/)

;Helper function
;a trial method for nearest neighbour interpolation in cases where you need an irregular grid
; Note: this is the same as the function nearestneighbor in tinterpol_mxn.pro. If you find a bug here
; please fix the other version too.
function spd_ui_nearestneighbor, v, x, u
; v: these are the actual values that will be interpolated (interpolates along one dimension)
; x: these are the x values (probably time) corresponding to the data (v) values (one dim array)
; u: this is the new array of x values, again 1 dim

n_u = n_elements(u)
n_x = n_elements(x)

;-------------------------------------------------------------------------------------
;Method 1: should work even if x is not monotonically incr. But can be very slow or
;run out of memory for large arrays
;-------------------------------------------------------------------------------------
;;this puts into the variable index the index of the closest value in x to each value in u
;; the resulting index is a one-dim subscript of a 2-dim array and thus must be manipulated further
;mindiff = min(abs(rebin(x,n_x,n_u) - rebin(transpose(u),n_x,n_u)),index,dimension=1)
;; convert index to multidim subscript
;index2d = array_indices([n_x,n_u],index,/dimensions)
;; only the first column of index2d is useful, containing index into x of nearest neigbor to each point in u
;; second column is just 0,1,2,3,..etc
;actualindex = transpose(index2d[0,*])

;-------------------------------------------------------------------------------------
;Method 2: only works if x is monotonically increasing or decreasing (same is true for interpol)
;Should be faster than above method for large arrays.
;-------------------------------------------------------------------------------------
;value_locate brackets each u in x
nearvalue = transpose(value_locate(x,u))
neararray = [nearvalue>0, (nearvalue+1)<(n_x-1)]; form an array with columns nearvalue and nearvalue +1, but restrict to valid indices into n_x
mindiff = min(abs(x[neararray]-rebin(transpose(u),2,n_u)),index, dimension=1)
; index gives you the (1 dim) index in neararray, indicating whether nearvalue or nearvalue+1 is closer
; the mod converts simply to 0 or 1
actualindex = transpose(nearvalue>0 + (index mod 2))


output = v[actualindex]
return, output
end

;helper function
;Does the actual interpolation on one N-dimensional quantity(but only interpolating along dimension 1)
;Function is originally from tinterpol_mxn
function spd_ui_interpolate_vec_mxn,v,x,u,nearest_neighbor=nearest_neighbor,_extra=_extra

compile_opt idl2, hidden

n = n_elements(u)

if n le 0 then return,-1L

;if the value is atomic return it
if(size(v,/n_dim) eq 0) then begin 
    error=1
    return,replicate(v,n_elements(u))
endif

v_s = size(v,/dimensions)

;handle single input case...it should extrapolate a constant matrix
if(v_s[0] eq 1) then begin
    v_s[0] = 2
    v = rebin(v,v_s)
    x = rebin([x],2)
    x[1] = x[0] + 1.0 ;so the timeseries ascends
endif

;I think I actually handled the 1 case generally
;if(n_elements(v_s) eq 1) then return,interpol(v,n)

v_s_o = v_s

v_s_o[0] = n

out = dindgen(v_s_o)

;the transpose and the reverse make the indexing scheme work out
;cause the in variables(and tplot variables) work more or less in row
;row major, but idl indexes column major
out_idx = transpose(lindgen(reverse(v_s_o)))

in_idx = transpose(lindgen(reverse(v_s)))

;calculate the number of elements in each matrix/vectors/whatever

product = 1

if n_elements(v_s) gt 1 then begin
  product = product(v_s[1:*])
endif

;for i = 1,n_elements(v_s) - 1L do begin
;
;    product *= v_s[i]
;
;endfor

for i = 0,product-1L do begin

    idx1 = where((out_idx mod product) eq i)
    idx2 = where((in_idx mod product) eq i)

    if(size(idx1,/n_dim) eq 0 || $
       n_elements(idx1) ne n || $
       size(idx2,/n_dim) eq 0 || $
       n_elements(idx2) ne v_s[0]) $
       then return, -1L

    if not keyword_set(u) then begin
      if keyword_set(nearest_neighbor) then begin
        out[idx1] = congrid( v[idx2],n)
      endif else out[idx1] = interpol(v[idx2],n,_extra=_extra) 
    endif else begin
      if keyword_set(nearest_neighbor) then begin
        out[idx1] = spd_ui_nearestneighbor(v[idx2],x,u)
      endif else out[idx1] = interpol(v[idx2],x,u,_extra=_extra)
    endelse

endfor
; for nearest neighbor case cast the output type to the input type
; This is so that if you interpolate bit-packed data you will still be
; able to use bitplot to plot it. It may be that all data should be cast
; to its input type - or maybe not.
if keyword_set(nneigbor) then begin
  out = fix(out, type=size(v, /type))
endif
return,out

end



;+
;NAME:
;  spd_ui_interpolate
;
;PURPOSE:
;  Interpolates over x-axis of active data and adds new data to
;  loaded data object.  Intended to be called from spd_ui_dproc.
;
;CALLING SEQUENCE:
;  spd_ui_interpolate, result, loadedData, historywin, statusbar
;
;INPUT:
;  loadedData: Loaded data object.
;  historywin: History window reference.
;  statusbar: Status bar reference.
;  guiid: If set and valid, will make user queries part of GUI widget hierarchy.
;  result: Anonymous structure returned by spd_ui_interpol_options:
;          {
;           num: Number of points for the result when not matching.
;           cad: Cadence of the result when using that option
;           match: Flag indicating to interpolate to another variable's abcissa.
;           matchto: Group name of variable to match to.
;           extrap: Flags for extrapolation type: [none, last value, NaN]
;           suffix: Suffix to create new data group with.
;           limit: Flag indicating time limits.
;           trange: Time range array containing desired time limits.
;           type: Flags for interpolation type: [linear,quadratic,lsquadratic,spline]
;           ntype: Flag to use number or cadence (0,1 respectively)
;           ok: (not used) Indicates spd_ui_interpol_options was successful
;           }
;
;HISTORY:
;
;-

pro spd_ui_interpolate, result,in_vars, loadedData, historywin, statusbar, $
                        fail=fail, guiid=guiid, _extra=_extra, replay=replay, $
                        overwrite_selections = overwrite_selections, $
                        cadence_selections = cadence_selections
    compile_opt idl2

  catch, on_err
  if on_err ne 0 then begin
    catch, /cancel
    help, /last_message, output=msg
    for i=0, n_elements(msg)-1 do historywin->update,msg[i]
    ok = error_message('Error in Interpolate function. Interpolation halted.', $
                       /center,/noname,title='Interpolation Error')
    fail=1
    return
  endif
  
  if ~keyword_set(guiId) then begin
    guiId = 0
  endif


active_data = in_vars


;Initializations
fail=0b
new_active=''
skipped=''
_extra = {QUADRATIC:result.type[1], $
          LSQUADRATIC:result.type[2], $
          SPLINE:result.type[3]}
nearest_neighbor = result.type[4] 

overwrite_selection = ''
cadence_selection= ''

overwrite_count = 0
cadence_count = 0

if undefined(replay) then begin
    overwrite_selections = ''
    cadence_selections = ''
endif

;Loop over active data
for j=0, n_elements(active_data)-1 do begin

  ;Get data elements
  loadedData->getvardata, name=active_data[j], $
                          time = t, $
                          data = d, $
                          yaxis = yd, $
                          dlimit=dl, $
                          limits=l
      
  if ~ptr_valid(t) || ~ptr_valid(d) then begin
    warning_message = active_data[j] + 'is invalid please re-try'
    ok = spd_ui_prompt_widget(guiId,statusbar,historywin,promptText=warning_message,title='Interpolate error', frame_attr=8)
    ;fail = 1
    continue
  endif 
      
  if ptr_valid(dl) then begin
    dlimit = *dl
  endif 
  
  if ptr_valid(l) then begin
    limit = *l
  endif

  ydata = ptr_valid(yd) ? 1b:0b

  ;Get time restrictions
  n_t = n_elements(*t)
  if result.limit[0] ne 0 then begin
    t0 = result.trange[0]
    t1 = result.trange[1]
  endif else begin
    t0 = (*t)[0]
    t1 = (*t)[n_t-1]
  endelse
  
  if result.limit[0] ne 0 then begin
    ;Apply time limits 
    idx = where( ((*t) ge result.trange[0]) and ((*t) le result.trange[1]), c)
    if c eq 0 then begin
      warning_message = 'No data in selected range for ' + active_data[j] + ' please re-try'
      ok = spd_ui_prompt_widget(guiId,statusbar,historywin,promptText=warning_message,title='Interpolate error', frame_attr=8)
      ;fail = 1
      continue
    endif
    ;Ensure first point ousite the time range is included
    if idx[0] gt 0 then idx = [ idx[0]-1, [idx] ]
    if idx[n_elements(idx)-1] lt (n_t-1) then idx = [ [idx], idx[n_elements(idx)-1]+1 ]  
  endif else begin
    idx = lindgen(n_t)
  endelse


  ;Check for matching
  ;----------
    if result.match[0] ne 0 then begin

;      if result.matchto[0] eq active_data[j] then begin
;        statusbar->update,'Skipping '+active_data[j]+', cannot match to same quantity.'
;        continue
;      endif else 
      
      statusbar->update,'Interpolate: Matching to '+result.matchto[0]+'...'

      ;Get new abscissa for interpol function
      loadedData->getvardata, name=result.matchto[0], time = u
     
      ;apply limit to matching variable if user limit is set, or if extrapolation is turned off
      if result.extrap[0] ne 0 || result.limit[0] ne 0 then begin
      
       ;Apply time limits to result's abcissa
        idxu = where( ((*u) ge t0) and ((*u) le t1), c)
        if c eq 0 then begin
          ok = error_message('No data to match in selected range, please re-try', $
                       /center,/noname,traceback=0, title='Interpolate error')
  ;        fail = 1
          skipped = [skipped,active_data[j]]
          continue
        endif
        
        t_p = (*u)[idxu]
      endif else begin
        t_p = (*u)
      endelse
      
      d_p = spd_ui_interpolate_vec_mxn(*d,*t,t_p,nearest_neighbor=nearest_neighbor,_extra=_extra)
      
      if ydata eq 1 then begin
       yd_p = spd_ui_interpolate_vec_mxn(*yd,*t,t_p,nearest_neighbor=nearest_neighbor,_extra=_extra)
      endif
        
      if result.extrap[0] eq 0 then begin
      
        ;find out of range indices to be extrapolated
        if result.limit[0] ne 0 then begin     
          low_idx = where(t_p lt t0,low_c)
          high_idx = where(t_p gt t1,high_c)
        endif else begin
          low_idx = where(t_p lt (*t)[0],low_c)
          high_idx =  where(t_p gt (*t)[n_t-1],high_c)
        endelse
        
        ;extrapolate low end indices
        if low_c gt 0 then begin
          if result.extrap[1] ne 0 then begin ;extrapolate by repetition
            d_p[low_idx,*] = (*d)[0,*] ## (dblarr(low_c)+1)
            if ydata eq 1 then begin
              yd_p[low_idx,*] = (*yd)[0,*] ## (dblarr(low_c)+1)
            endif
          endif else if result.extrap[2] ne 0 then begin ;extrapolate with nans
            d_p[low_idx,*] = !VALUES.D_NAN
;            if ydata eq 1 then begin ;we probably don't want this
;              yd_p[low_idx,*] = !VALUES.D_NAN
;            endif
          endif
        endif
        
        ;extrapolate high end indices
        if high_c gt 0 then begin
          if result.extrap[1] ne 0 then begin ;extrapolate by repetition
            d_p[high_idx,*] = (*d)[n_t-1,*] ## (dblarr(high_c)+1)
            if ydata eq 1 then begin
              yd_p[high_idx,*] = (*yd)[n_t-1,*] ## (dblarr(high_c)+1)
            endif
          endif else if result.extrap[2] ne 0 then begin ;extrapolate with nans
            d_p[high_idx,*] = !VALUES.D_NAN
;            if ydata eq 1 then begin ;we probably don't want this
;              yd_p[high_idx,*] = !VALUES.D_NAN
;            endif
          endif
        endif
        
      endif
      
  
  ;Check if using Cadence
  ;---------
    endif else if result.ntype[0] then begin


      statusbar->update,'Interpolate: Using '+strtrim(result.cad[0])+' second cadence...'


      n_d0 = n_elements((*d)[idx,0])
      n_d1 = n_elements((*d)[0,*])
      if ydata then begin
        n_yd1 = n_elements((*yd)[0,*])
      endif

      ;Create new abcissa at desired cadence
      n_u = floor( (t1-t0)/result.cad[0],/l64)

      u = result.cad[0]*dindgen(n_u) + t0

      ;Check requested resolution
      if n_d0 gt n_u then begin
          if cadence_selection ne 'yestoall' then begin
                if cadence_selection eq 'notoall' then continue
                
                cadence_selection = ''
                if ~undefined(replay) then begin
                    if cadence_count gt n_elements(cadence_selections) then begin
                        dprint, dlevel = 0, 'Error: discrepancy in SPEDAS document, may have led to a document load error'
                        cadence_selection = 'yestoall'
                    endif else begin
                        cadence_selection = cadence_selections[cadence_count]
                    endelse
                endif else begin
                    prompt = 'Using the requested cadence  ('+strtrim(result.cad[0],2)+ $
                             ' sec) will decrease time resolution for '+active_data[j]+ $
                             ' .  Proceed?'
                    cadence_selection = spd_ui_prompt_widget(guiId,statusbar,historywin,promptText=prompt,$
                                              /yes,/no,/allyes,/allno,title='Decrease time resolution?', frame_attr=8)
                    cadence_selections = array_concat_wrapper(cadence_selection, cadence_selections)
                endelse
                
                cadence_count++
                if (cadence_selection eq 'no' || cadence_selection eq 'notoall') then continue
          endif
      endif

      ;Do interpolation
      ;-----------
      d_p = dblarr( n_u, n_d1, /nozero)
      for i=0, n_d1-1 do begin
        if nearest_neighbor eq 1 then begin 
          d_p[0,i] = congrid( (*d)[idx,i],n_u)
        endif else d_p[0,i] = interpol( (*d)[idx,i], (*t)[idx], (u), _extra=_extra )
      endfor

      ;Interpolate over any yaxis data
      ;-----------
      if ydata then begin
        yd_p = dblarr( n_u, n_yd1, /nozero)
        for i=0, n_yd1-1 do begin
          if nearest_neighbor eq 1 then begin
            yd_p[0,i] = congrid( (*yd)[idx,i],n_u)
          endif else yd_p[0,i] = interpol( (*yd)[idx,i], (*t)[idx], (u), _extra=_extra )
        endfor
      endif

      t_p = temporary(u)

  ;If not, use default option
  ;---------
    endif else begin

      statusbar->update,'Interpolate: Using '+strtrim(result.num[0])+' points...'

      n_d0 = n_elements((*d)[idx,0])
      n_d1 = n_elements((*d)[0,*])
      if ydata then begin
        n_yd1 = n_elements((*yd)[0,*])
      endif

      ;Create new x-axis explicitly
      u = interpol( [t0,t1], result.num[0])
      
      ;Check requested resolution
      if n_d0 gt result.num[0] then begin
          if cadence_selection ne 'yestoall' then begin
                if cadence_selection eq 'notoall' then continue
                
                cadence_selection = ''
                if ~undefined(replay) then begin
                    if cadence_count gt n_elements(cadence_selections) then begin
                        dprint, dlevel = 0, 'Error: discrepancy in SPEDAS document, may have led to a document load error'
                        cadence_selection = 'yestoall'
                    endif else begin
                        cadence_selection = cadence_selections[cadence_count]
                    endelse
                endif else begin
                    prompt = 'The requested number of data points is less than '+$
                             'the current data resolution for '+active_data[j]+ $
                             ' ( '+strtrim(n_d0,2)+' points).  Proceed?'
                    cadence_selection = spd_ui_prompt_widget(guiId,statusbar,historywin,promptText=prompt,$
                                              /yes,/no,/allyes,/allno,title='Decrease time resolution?', frame_attr=8)
                    cadence_selections = array_concat_wrapper(cadence_selection, cadence_selections)
                endelse
                
                cadence_count++
                if (cadence_selection eq 'no' || cadence_selection eq 'notoall') then continue
          endif
      endif

      ;Do interpolation
      ;-----------
      d_p = dblarr(result.num[0], n_d1, /nozero)
      for i=0, n_d1-1 do begin
        if nearest_neighbor eq 1 then begin
          d_p[0,i] = congrid( (*d)[idx,i], result.num[0])
        endif else d_p[0,i] = interpol( (*d)[idx,i], (*t)[idx], u, _extra=_extra)
      endfor

      ;Interpolate over any yaxis data
      ;----------- 
      if ydata then begin
        yd_p = dblarr(result.num[0], n_yd1, /nozero)
        for i=0, n_yd1-1 do begin
          if nearest_neighbor eq 1 then begin
            yd_p[0,i] = congrid( (*yd)[idx,i], result.num[0])
          endif else yd_p[0,i] = interpol( (*yd)[idx,i], (*t)[idx], u, _extra=_extra)
        endfor
      endif

      t_p = temporary(u)

    endelse


    ;Set up data to be added
    name = active_data[j]+result.suffix[0]
    if ptr_valid(l) then limit = *l else limit=0 
    if ptr_valid(dl) then dlimit = *dl else dlimit=0
    if ydata then data = {x:temporary(t_p), y:temporary(d_p), v:temporary(yd_p)} $
      else data = {x:temporary(t_p), y:temporary(d_p)}

    ; check if the new tplot variable already exists, query the user to overwrite it if it does
    spd_ui_check_overwrite_data,name,loadedData,guiId,statusBar,historyWin,overwrite_selection,overwrite_count,$
                             replay=replay,overwrite_selections=overwrite_selections
    if strmid(overwrite_selection, 0, 2) eq 'no' then continue

    ; Add data
    add = loadedData->addData(name, data, dlimit=dlimit, limit=limit, isspect=ydata)
    if add then new_active = [new_active,name] $
      else skipped = [skipped,active_data[j]]

endfor


;Reset active data quantities
  if n_elements(new_active) gt 1 then begin
    loadedData->clearAllActive
    for i=1, n_elements(new_active)-1 do loadedData->setactive,new_active[i]
  endif

;Return messages for any problems encountered
  if fail then begin
    statusbar->update,'Interpolate:  One or more problems occured. See history.'
  endif else begin
    statusbar->update,'Interpolate: Successful'
  endelse

  if n_elements(skipped) gt 1 then begin
    for i=1, n_elements(skipped)-1 do begin
      historywin -> update, skipped[i]+' not processed.'
    endfor
    statusbar->update, 'Interpolate: Some quantities were not processed. See history.'
  endif


end