[frames] | no frames]

# 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
<!--
expandto(location.href);
// -->

```

 Generated by Epydoc 3.0.1 on Thu Sep 22 04:25:12 2011 http://epydoc.sourceforge.net