Source Code for Module reflectometry.model1d.model.bspline

  1  # This program is public domain 
  2  """ 
  3  BSpline functions.  Given a set of knots, compute the 
  4  degree 3 Bspline and any derivatives that are required. 
  5  """ 
  6  import numpy 
  7   
  8  #=========================================================== 
  9  # BSpline Class 
 10  #=========================================================== 
 11 -class BSpline: 
 12      """ 
 13      Deprecated.  This class will be removed. 
 14      """ 
 15 -    def __init__(self, valList): 
 16          self._valList  = valList 
 17   
 18   
 19 -    def build(self): 
 20          return "BSpline(%s)" %( self._valList ) 
 21   
 22   
 23 -    def run(self): 
 24          return [self._valList ] 
 25   
 26   
 27 -    def runList(self): 
 28          ret = [] 
 29          for i in xrange( len(self._valList) ): 
 30              vals = "%.6f"%(self._valList[i]) 
 31              ret.append( vals ) 
 32   
 33          return ret 
 34   
 35 -    def run3(self): 
 36          lenList = len(self._valList) 
 37   
 38          ret = "[" 
 39          for i in xrange( lenList ): 
 40              ret += "%.3f"%(self._valList[i]) 
 41              if i != lenList-1: 
 42                  ret += "," 
 43          ret += "]" 
 44   
 45          return ret 
 46   
 47   
 48  bspline =  BSpline 
 49  spline  =  BSpline 
 50  slope   =  BSpline 
 51  Slope   =  BSpline 
 52   
 53   
 54   
 55 -def max(a,b): 
 56      return (a<b).choose(a,b) 
 57   
 58 -def min(a,b): 
 59      return (a>b).choose(a,b) 
 60   
 61   
 62 -def bspline3(knot, 
 63               control, 
 64               t, 
 65               clamp=True, 
 66               nderiv=0 
 67               ): 
 68      """ 
 69      Evaluate the B-spline specified by the given knot sequence and 
 70      control values at the parametric points t.   The knot sequence 
 71      should be four elements longer than the control sequence. 
 72   
 73      If nderiv is greater than 0, returns derivatives as well as the 
 74      function value.  For example, nderiv=3 returns f,f',f'',f''', 
 75   
 76      If clamp is True, the spline value is clamped to the value of 
 77      the final control point beyond the ends of the knot sequence. 
 78      If clamp is False, the spline will go to zero at +/- infinity 
 79      as in the traditional algorithm. 
 80      """ 
 81      degree = len(knot) - len(control); 
 82      if degree != 4: 
 83          raise ValueError, "must have two extra knots at each end" 
 84   
 85      if clamp: 
 86          # Alternative approach spline is clamped to initial/final control values 
 87          control = numpy.concatenate(([control[0]]*(degree-1), 
 88                                       control, [control[-1]])) 
 89      else: 
 90          # Traditional approach: spline goes to zero at +/- infinity. 
 91          control = numpy.concatenate(([0]*(degree-1), control, [0])) 
 92   
 93   
 94      # Deal with values outside the range 
 95      valid = (t > knot[0]) & (t <= knot[-1]) 
 96      tv  = t[valid] 
 97      f   = numpy.zeros(t.shape) 
 98      f[t<=knot[0]]  = control[0] 
 99      f[t>=knot[-1]] = control[-1] 
100   
101      # Find B-Spline parameters for the individual segments 
102      end     = len(knot)-1 
103      segment = knot.searchsorted(tv)-1 
104      tm2 = knot[max(segment-2,0)] 
105      tm1 = knot[max(segment-1,0)] 
106      tm0 = knot[max(segment-0,0)] 
107      tp1 = knot[min(segment+1,end)] 
108      tp2 = knot[min(segment+2,end)] 
109      tp3 = knot[min(segment+3,end)] 
110   
111      P4 = control[min(segment+3,end)] 
112      P3 = control[min(segment+2,end)] 
113      P2 = control[min(segment+1,end)] 
114      P1 = control[min(segment+0,end)] 
115   
116      # Compute second and third derivatives. 
117      if nderiv > 1: 
118          # Normally we require a recursion for Q, R and S to compute 
119          # df, d2f and d3f respectively, however Q can be computed directly 
120          # from intermediate values of P, S has a recursion of depth 0, 
121          # which leaves only the R recursion of depth 1 in the calculation 
122          # below. 
123          Q4 = (P4 - P3) * 3 / (tp3-tm0) 
124          Q3 = (P3 - P2) * 3 / (tp2-tm1) 
125          Q2 = (P2 - P1) * 3 / (tp1-tm2) 
126          R4 = (Q4 - Q3) * 2 / (tp2-tm0) 
127          R3 = (Q3 - Q2) * 2 / (tp1-tm1) 
128          if nderiv > 2: 
129              S4 = (R4 - R3) / (tp1-tm0) 
130              d3f = numpy.zeros(t.shape) 
131              d3f[valid] = S4 
132          R4 = ( (tv-tm0)*R4 + (tp1-tv)*R3 ) / (tp1 - tm0) 
133          d2f = numpy.zeros(t.shape) 
134          d2f[valid] = R4 
135   
136      # Compute function value and first derivative 
137      P4 = ( (tv-tm0)*P4 + (tp3-tv)*P3 ) / (tp3 - tm0) 
138      P3 = ( (tv-tm1)*P3 + (tp2-tv)*P2 ) / (tp2 - tm1) 
139      P2 = ( (tv-tm2)*P2 + (tp1-tv)*P1 ) / (tp1 - tm2) 
140      P4 = ( (tv-tm0)*P4 + (tp2-tv)*P3 ) / (tp2 - tm0) 
141      P3 = ( (tv-tm1)*P3 + (tp1-tv)*P2 ) / (tp1 - tm1) 
142      if  nderiv >= 1: 
143          df = numpy.zeros(t.shape) 
144          df[valid] = (P4-P3) * 3 / (tp1-tm0) 
145      P4 = ( (tv-tm0)*P4 + (tp1-tv)*P3 ) / (tp1 - tm0) 
146      f[valid]  = P4 
147   
148      if   nderiv == 0: return f 
149      elif nderiv == 1: return f,df 
150      elif nderiv == 2: return f,df,d2f 
151      else:             return f,df,d2f,d3f 
152