Source Code for Module ComboCode.cc.ivs.sigproc.interpol

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