Source Code for Module ivs.sigproc.interpol

  1  """ 
  2  Non-standard interpolation methods. 
  3  """ 
  4  import numpy as np 
  5  from scipy import ndimage 
  6  import astropy.io.fits as pf 
  7  import time 
  8  import itertools 
  9  import pyfinterpol 
 10   
 11 -def __df_dx(oldx,oldy,index,sharp=False): 
 12      """ 
 13      Preliminary estimation of df/dx 
 14      """ 
 15      xm1,fm1 = oldx[index-1],oldy[index-1] 
 16      x  ,f   = oldx[index]  ,oldy[index] 
 17      xp1,fp1 = oldx[index+1],oldy[index+1] 
 18       
 19      if not sharp: 
 20          df_dx = 1./(xp1-xm1) * ( fp1*(x-xm1)/(xp1-x) - fm1*(xp1-x)/(x-xm1)) + \ 
 21              f*(xp1-2*x+xm1) / ( (x-xm1)*(xp1-x)) 
 22      else: 
 23          df_dx = min( (fp1-f)/(xp1-x), (f-fm1)/(x-xm1)  ) 
 24      return df_dx 
 25   
 26  #-- interpolation by polynomial of degree 3 
 27 -def __P1(x,x0,x1): return  (x-x1)**2 * (2*x-3*x0+x1) / (x1-x0)**3 
 28 -def __P2(x,x0,x1): return -(x-x0)**2 * (2*x-3*x1+x0) / (x1-x0)**3 
 29 -def __P3(x,x0,x1): return  (x-x0)    * (x-x1)**2     / (x1-x0)**2 
 30 -def __P4(x,x0,x1): return  (x-x0)**2 * (x-x1)        / (x1-x0)**2 
 31   
 32 -def local_interpolation(newx,oldx,oldy,full_output=False): 
 33      """ 
 34      A local interpolation method by a polynomial of degree 3. 
 35       
 36      After Marc-Antoine Dupret, 2002 (extended version of PhD thesis). 
 37       
 38      Cannot extrapolate! 
 39       
 40      >>> np.random.seed(1114) 
 41      >>> oldx = np.sort(np.random.uniform(size=10)) 
 42      >>> oldy = oldx**2#np.random.uniform(size=10) 
 43      >>> oldy[5:] = -oldx[5:]**2 
 44      >>> newx = np.linspace(oldx.min(),oldx.max(),1000) 
 45      >>> newy,disconts = local_interpolation(newx,oldx,oldy,full_output=True) 
 46      >>> newy_ = local_interpolation(newx,oldx,oldy) 
 47       
 48      >>> sharpy = newy.copy() 
 49      >>> sharpy[disconts] = np.nan 
 50      >>> smoothy = newy.copy() 
 51      >>> smoothy[-disconts] = np.nan 
 52      >>> p = pl.figure() 
 53      >>> p = pl.plot(oldx,oldy,'ks-',ms=10) 
 54      >>> p = pl.plot(newx,smoothy,'go-',lw=2,ms=2,mec='g') 
 55      >>> p = pl.plot(newx,sharpy,'ro-',lw=2,ms=2,mec='r') 
 56      >>> p = pl.plot(newx,newy_,'bo--',lw=2,ms=2,mec='b') 
 57       
 58      @param newx: new x-array to interpolate on 
 59      @type newx: ndarray 
 60      @param oldx: old x-array to interpolate from 
 61      @type oldx: ndarray 
 62      @param oldy: old y-array to interpolate from 
 63      @type oldy: ndarray 
 64      @param full_output: also return points where discontinuities were detected 
 65      @type full_output: boolean 
 66      @return: interpolated array(, discontinuities) 
 67      @rtype: ndarray(,ndarray) 
 68      """ 
 69      #if not full_output: 
 70      #    return pyfinterpol.local_interpolation(newx,oldx,oldy) 
 71      #-- extend axis to be able to interpolate last point 
 72      lastx = 2*oldx[-1]-oldx[-2] 
 73      lasty = (oldy[-1]-oldy[-2])/(oldx[-1]-oldx[-2])*(oldx[-1]-lastx) + oldy[-1] 
 74      oldy = np.hstack([oldy,lasty]) 
 75      oldx = np.hstack([oldx,lastx]) 
 76      #-- prepare new y array 
 77      newy = np.zeros_like(newx) 
 78      #-- keep information on where sharp features occur, if asked 
 79      if full_output: 
 80          disconts = np.zeros(len(newx),bool) 
 81          #disconts = np.zeros(len(newx)) 
 82       
 83      index = -1 
 84      for i,x in enumerate(newx): 
 85          index = oldx.searchsorted(x) 
 86           
 87          #if index>=(len(oldx)-1): continue 
 88          x0,f0 = oldx[index-1],oldy[index-1] 
 89          x1,f1 = oldx[index],oldy[index] 
 90          x2,f2 = oldx[index-2],oldy[index-2] 
 91           
 92          #-- check sharpness of feature     
 93          #sharpness_ = 1./((f1-f0)/(x1-x0)*(x1-x2)/(f1-f2)) 
 94          numerator = ((x1-x0)*(f1-f2)) 
 95          denominator = ((f1-f0)*(x1-x2)) 
 96          #print 'p',numerator,denominator 
 97          if denominator==0: 
 98              sharpness = 0. 
 99          else: 
100              sharpness = numerator/denominator 
101          sharp = (0.2<=sharpness<=0.5)  
102          if full_output: 
103              disconts[i] = sharp#ness#sharp 
104          #-- preliminary estimation of df/dx 
105          dfdx0 = __df_dx(oldx,oldy,index-1,sharp=sharp) 
106          dfdx1 = __df_dx(oldx,oldy,index,sharp=sharp) 
107           
108           
109          #-- interpolation by polynomial of degree 3 
110          P1 =  (x-x1)**2 * (2*x-3*x0+x1) / (x1-x0)**3 
111          P2 = -(x-x0)**2 * (2*x-3*x1+x0) / (x1-x0)**3 
112          P3 =  (x-x0)    * (x-x1)**2     / (x1-x0)**2 
113          P4 =  (x-x0)**2 * (x-x1)        / (x1-x0)**2 
114           
115          #-- interpolating polynomial 
116          Px = f0*P1 + f1*P2 + dfdx0*P3 + dfdx1*P4 
117          newy[i] = Px 
118      if full_output: 
119          return newy,disconts 
120      else: 
121          return newy 
122   
123   
124 -def local_interpolation_ND(newx,oldx,oldy): 
125      #-- oldx should be a list of [x,y,z,...] 
126      #   oldy should be array of N(x) x N(y) x N(z) x ... 
127      points = np.zeros((len(args),3)) 
128      for _x in args: 
129          index = oldx.searchsorted(_x) 
130          x0,f0 = oldx[index-1],oldy[index-1] 
131          x1,f1 = oldx[index],oldy[index] 
132          x2,f2 = oldx[index-2],oldy[index-2] 
133       
134      polynomials = [] 
135       
136   
137 -def create_pixeltypegrid(grid_pars,grid_data): 
138      """ 
139      Creates pixelgrid and arrays of axis values. 
140       
141      Starting from: 
142          * grid_pars: 2D numpy array, 1 column per parameter, unlimited number of cols 
143          * grid_data: 2D numpy array, 1 column per variable, data corresponding to the rows in grid_pars 
144       
145      The grid should be rectangular and complete, i.e. every combination of the unique values in the  
146      parameter columns should exist. If not, a nan value will be inserted. 
147       
148      @param grid_pars: Npar x Ngrid array of parameters 
149      @type grid_pars: array 
150      @param grid_data: Ndata x Ngrid array of data 
151      @type grid_data:array 
152      @return: axis values and pixelgrid 
153      @rtype: array, array 
154      """ 
155   
156      uniques = [np.unique(column, return_inverse=True) for column in grid_pars] 
157      #[0] are the unique values, [1] the indices for these to recreate the original array 
158   
159      axis_values = [uniques_[0] for uniques_ in uniques] 
160      unique_val_indices = [uniques_[1] for uniques_ in uniques] 
161       
162      data_dim = np.shape(grid_data)[0] 
163   
164      par_dims   = [len(uv[0]) for uv in uniques] 
165   
166      par_dims.append(data_dim) 
167      pixelgrid = np.ones(par_dims) 
168       
169      # We put np.inf as default value. If we get an inf, that means we tried to access 
170      # a region of the pixelgrid that is not populated by the data table 
171      pixelgrid[pixelgrid==1] = np.inf 
172       
173      # now populate the multiDgrid 
174      indices = [uv[1] for uv in uniques] 
175      pixelgrid[indices] = grid_data.T 
176      return axis_values, pixelgrid 
177   
178 -def interpolate(p, axis_values, pixelgrid): 
179      """ 
180      Interpolates in a grid prepared by create_pixeltypegrid(). 
181       
182      p is an array of parameter arrays 
183       
184      @param p: Npar x Ninterpolate array 
185      @type p: array 
186      @return: Ndata x Ninterpolate array 
187      @rtype: array 
188      """ 
189      # convert requested parameter combination into a coordinate 
190      #p_ = [np.searchsorted(av_,val) for av_, val in zip(axis_values,p)] 
191      # we force the values to be inside the grid, to avoid edge-effect rounding 
192      # (e.g. 3.099999 is edge, while actually it is 3.1). For values below the 
193      # lowest value, this is automatically done via searchsorted (it return 0) 
194      # for values higher up, we need to force it 
195      p_ = [] 
196      for av_,val in zip(axis_values,p): 
197          indices = np.searchsorted(av_,val) 
198          indices[indices==len(av_)] = len(av_)-1 
199          p_.append(indices) 
200       
201   
202      #-- The type of p is changes to the same type as in axis_values to catch possible rounding errors 
203      #   when comparing float64 to float32. 
204      for i, ax in enumerate(axis_values): 
205          p[i] = np.array(p[i], dtype = ax.dtype) 
206       
207      #-- Convert requested parameter combination into a coordinate 
208      p_ = np.array([np.searchsorted(av_,val) for av_, val in zip(axis_values,p)]) 
209      lowervals_stepsize = np.array([[av_[p__-1], av_[p__]-av_[p__-1]] \ 
210                              for av_, p__ in zip(axis_values,p_)]) 
211      p_coord = (p-lowervals_stepsize[:,0])/lowervals_stepsize[:,1] + np.array(p_)-1 
212   
213      # interpolate 
214      return np.array([ndimage.map_coordinates(pixelgrid[...,i],p_coord, order=1, prefilter=False) \ 
215                  for i in range(np.shape(pixelgrid)[-1])]) 
216   
217   
218   
219   
220  if __name__=='__main__': 
221      import pylab as pl 
222      #import time 
223      #from doctest import testmod 
224      #testmod() 
225       
226      np.random.seed(1114) 
227      oldx = np.sort(np.random.uniform(size=10)) 
228      oldy = oldx**2#np.random.uniform(size=10) 
229      oldy[5:] = -oldx[5:]**2 
230      newx = np.linspace(oldx.min(),oldx.max(),1000) 
231      newy,disconts = local_interpolation(newx,oldx,oldy,full_output=True) 
232      newy_ = local_interpolation(newx,oldx,oldy) 
233       
234      sharpy = newy.copy() 
235      sharpy[disconts] = np.nan 
236      smoothy = newy.copy() 
237      smoothy[-disconts] = np.nan 
238      p = pl.figure() 
239      p = pl.plot(oldx,oldy,'ks-',ms=10) 
240      p = pl.plot(newx,smoothy,'go-',lw=2,ms=2,mec='g') 
241      p = pl.plot(newx,sharpy,'ro-',lw=2,ms=2,mec='r') 
242      p = pl.plot(newx,newy_,'bo--',lw=2,ms=2,mec='b') 
243       
244      pl.show() 
245