Source Code for Module gmisclib.g_lin_fit

  1  """Generic linear best fits.""" 
  2   
  3  import types 
  4  from gmisclib import Num 
  5   
  6   
  7 -def fit(idata, fitting_fcns, ff_info=None, wt=None): 
  8          """Fits  data to a linear combination of fitting_fcns. 
  9          @rtype: tuple 
 10          @return: The returned value is a tuple of: 
 11                  the coefficients of the optimal solution (as an array), 
 12                  the residual (a single float number), 
 13                  the rank of the fit (int), 
 14                  an array of the singular values, 
 15                  and the best fit to the data. 
 16                  The returned value follows Num.LA.linear_least_squares(), 
 17                  except that residual is always calculated and is a scalar. 
 18          @type idata: an array of (vector, info). 
 19          @type fitting_fcns: an array of functions. 
 20          @param fitting_fcns: the basis functions for the linear regression. 
 21                  The fitting functions return vectors (which must match the data). 
 22                  Fitting functions are called as f(info, ff_info, len(vector)). 
 23                  Vectors may have different lengths, as long as the 
 24                  corresponding data and fit vectors agree. 
 25          """ 
 26   
 27          nd = len(idata) 
 28          assert wt is None or len(wt) == nd 
 29          w = [None] * nd 
 30   
 31          # We concatenate the data vectors to do one big fit all at once. 
 32          # d_idx keeps track of where each vector starts and ends. 
 33          d_idx = Num.zeros((nd+1,), Num.Int) 
 34          nd_total = 0 
 35          for i in range(nd): 
 36                  assert len(idata[i]) == 2 
 37                  assert type(idata[i]) == types.TupleType 
 38                  d_idx[i] = nd_total             # Where does that dataset start? 
 39                  nd_total = nd_total + len(idata[i][0]) 
 40          d_idx[nd] = nd_total 
 41   
 42          # Build the data array: 
 43          data = Num.zeros((nd_total,), Num.Float) 
 44          for i in range(nd): 
 45                  if wt is None: 
 46                          w[i] = Num.ones((len(idata[i][0]),), Num.Float) 
 47                  else: 
 48                          assert len(wt[i]) == len(idata[i][0]) 
 49                          w[i] = Num.asarray(wt[i], Num.Float) 
 50                  assert len(idata[i][0]) == d_idx[i+1]-d_idx[i] 
 51                  data[d_idx[i] : d_idx[i+1]] = Num.asarray(idata[i][0], 
 52                                                                  Num.Float) * w[i] 
 53   
 54          # Build the array of basis functions: 
 55          nff = len(fitting_fcns) 
 56          base = Num.zeros((nd_total, nff), Num.Float) 
 57          for i in range(nd): 
 58                  for j in range(nff): 
 59                          ff = fitting_fcns[j] 
 60                          zfit = Num.asarray(ff(idata[i][1], ff_info, d_idx[i+1]-d_idx[i]), 
 61                                                  Num.Float) 
 62                          assert len(zfit.shape)==1 
 63                          assert zfit.shape[0] == d_idx[i+1]-d_idx[i] 
 64                          base[d_idx[i]:d_idx[i+1], j] = zfit * w[i] 
 65   
 66          rcond = 1e-5 
 67          x, resid, rank, s = Num.LA.linear_least_squares(base, data, rcond) 
 68          bfit = Num.matrixmultiply(base, x) 
 69          resid = Num.sum(Num.square(bfit-data).ravel()) 
 70   
 71          bestfit = [] 
 72          # Chop the fit array back up to match the input vectors: 
 73          for i in range(nd): 
 74                  bestfit.append(bfit[d_idx[i]:d_idx[i+1]]/w[i]) 
 75   
 76          return x, resid, rank, s, bestfit 
 77   
 78   
 79   
 80   
 81   
 82   
 83   
 84 -def test(): 
 85          test1() 
 86          test2() 
 87          test3() 
 88          test4() 
 89          test5() 
 90          test6() 
 91   
 92 -def test1(): 
 93          x, resid, rank, s, f = fit([([0,1,2,3,4,5], None)], [lambda di, fi, l: range(l)], None) 
 94          assert abs(x[0]-1)<0.0001 
 95          assert resid < 0.001 
 96          assert rank == 1 
 97   
 98 -def test2(): 
 99          x, resid, rank, s, f = fit([([0,1,2,3], None), ([0,1,2,3,4], None)], [lambda di, fi, l: range(l)], None) 
100          assert abs(x[0]-1)<0.0001 
101          assert resid < 0.001 
102          assert rank == 1 
103   
104 -def test3(): 
105          x, resid, rank, s, f = fit([([0,1,2,3], None), ([0,2,4,6], None)], [lambda di, fi, l: range(l)]) 
106          assert abs(x[0]-1.5)<0.0001 
107          assert rank == 1 
108          x, resid, rank, s, f = fit([([0,1,2,3], None), ([0,2,4,6], None)], [lambda di, fi, l: range(l)], wt=[[1,1,1,1],[1e-10,1e-10,1e-10,1e-10]]) 
109          assert abs(x[0]-1)<0.0001 
110          assert rank == 1 
111   
112 -def test4(): 
113          x, resid, rank, s, f = fit([([0,1,2,3], None), ([1,2,3,4], None)], [lambda di, fi, l: range(l), lambda di, fi, l: [1]*l], None) 
114          assert abs(x[0]-1)<0.0001 
115          assert abs(x[1]-0.5)<0.0001 
116          assert rank == 2 
117   
118 -def test5(): 
119          x, resid, rank, s, f = fit([([1,2,3,4], None)], [lambda di, fi, l: range(l), lambda di, fi, l: [1]*l]) 
120          assert abs(x[0]-1)<0.0001 
121          assert abs(x[1]-1)<0.0001 
122          assert resid < 0.001 
123          assert rank == 2 
124          assert abs(f[0][0]-1)<0.001 
125          assert abs(f[0][1]-2)<0.001 
126   
127 -def test6(): 
128          x, resid, rank, s, f = fit( 
129                                          [([1,2,3,4], None)], 
130                                          [lambda di, fi, l: range(l), lambda di, fi, l: [1]*l], 
131                                          wt=[[5,2,3,1e-10]] 
132                                  ) 
133          assert abs(x[0]-1)<0.0001 
134          assert abs(x[1]-1)<0.0001 
135          assert resid < 0.001 
136          assert rank == 2 
137   
138   
139   
140   
141   
142  if __name__ == '__main__' : 
143          test() 
144          print "OK: passed tests" 
145