|
1 #!/usr/bin/python |
|
2 # Contains all the polynomial functions |
|
3 |
|
4 import pylab as pl |
|
5 |
|
6 def polsize(Q, degQ): |
|
7 """ Determines dimensions of a polynomial matrix. """ |
|
8 rQ, cQ = pl.atleast_2d(Q).shape |
|
9 cQ = cQ/float(degQ+1) |
|
10 if abs(round(cQ)-cQ) > 1e-6: |
|
11 print "Degree of input inconsistent with number of columns" |
|
12 return |
|
13 else: |
|
14 cQ = int(round(cQ)) |
|
15 return rQ, cQ |
|
16 |
|
17 def polmul(A, degA, B, degB): |
|
18 A = pl.atleast_2d(A) |
|
19 B = pl.atleast_2d(B) |
|
20 rA, cA = polsize(A, degA) |
|
21 rB, cB = polsize(B, degB) |
|
22 if cA != rB: |
|
23 print "polmul: Inconsistent dimensions of input matrices" |
|
24 return |
|
25 degC = degA + degB |
|
26 C = [] |
|
27 for k in range(0, degC+1): |
|
28 mi = 0 |
|
29 if k-degB > mi: |
|
30 mi = k-degB |
|
31 ma = degA |
|
32 if k < ma: |
|
33 ma = k |
|
34 Ck = pl.zeros((rA,cB)) |
|
35 for i in range(mi, ma+1): |
|
36 Ck = Ck + pl.dot(A[..., i*cA:(i+1)*cA], B[..., (k-i)*cB:(k-i+1)*cB]) |
|
37 Ck = pl.squeeze(Ck) |
|
38 C = pl.hstack((C, Ck)) |
|
39 return C, degC |
|
40 |
|
41 def poladd(A, degA, B, degB): |
|
42 A = pl.atleast_2d(A) |
|
43 B = pl.atleast_2d(B) |
|
44 rA, cA = polsize(A, degA) |
|
45 rB, cB = polsize(B, degB) |
|
46 if cA != rB: |
|
47 print "polmul: Inconsistent dimensions of input matrices" |
|
48 return |
|
49 degC = max(degA, degB) |
|
50 |
|
51 if degC >= degA: |
|
52 A = pl.hstack((A, pl.zeros((rA,(degC-degA)*cA)))) |
|
53 |
|
54 if degC >= degB: |
|
55 B = pl.hstack((B, pl.zeros((rB,(degC-degB)*cB)))) |
|
56 |
|
57 C = A + B |
|
58 return C, degC |
|
59 |
|
60 def polsplit2(fac, a=1-1e-5): |
|
61 fac = pl.atleast_1d(fac) |
|
62 if a>1: |
|
63 print "good polynomial also is unstable" |
|
64 return |
|
65 roots = pl.roots(fac) |
|
66 |
|
67 # extract good and bad roots |
|
68 badindex = pl.find(pl.absolute(roots)>=a-1.0e-5) |
|
69 badpoly = pl.poly(roots[badindex]) |
|
70 goodindex = pl.find(pl.absolute(roots)<a-1.0e-5) |
|
71 goodpoly = pl.poly(roots[goodindex]) |
|
72 # scale by equating the largest terms |
|
73 index = pl.absolute(fac).argmax() |
|
74 goodbad = pl.convolve(goodpoly, badpoly) |
|
75 factor = fac[index]/goodbad[index] |
|
76 goodpoly = goodpoly * factor |
|
77 badpoly = pl.atleast_1d(badpoly) |
|
78 goodpoly = pl.atleast_1d(goodpoly) |
|
79 return goodpoly, badpoly |
|
80 |
|
81 def polsplit3(fac, a=1): |
|
82 fac = pl.atleast_1d(fac) |
|
83 if a>1: |
|
84 print "good polynomial also is unstable" |
|
85 return |
|
86 roots = pl.roots(fac) |
|
87 |
|
88 # extract good and bad roots |
|
89 badindex = pl.find((pl.absolute(roots)>=a-1.0e-5) + (pl.real(roots)<-0.05)) |
|
90 badpoly = pl.poly(roots[badindex]) |
|
91 goodindex = pl.find((pl.absolute(roots)<a-1.0e-5) * (pl.real(roots)>=-0.05)) |
|
92 goodpoly = pl.poly(roots[goodindex]) |
|
93 # scale by equating the largest terms |
|
94 index = pl.absolute(fac).argmax() |
|
95 goodbad = pl.convolve(goodpoly, badpoly) |
|
96 factor = fac[index]/goodbad[index] |
|
97 goodpoly = goodpoly * factor |
|
98 badpoly = pl.atleast_1d(badpoly) |
|
99 goodpoly = pl.atleast_1d(goodpoly) |
|
100 return goodpoly, badpoly |
|
101 |
|
102 def putin(A, degA, B, degB, i, j): |
|
103 from clcoef import clcoef |
|
104 A = pl.atleast_2d(A) |
|
105 B = pl.atleast_2d(B) |
|
106 rA, cA = polsize(A,degA) |
|
107 if degB > degA: |
|
108 A = pl.hstack((A, pl.zeros((rA,(degB-degA)*cA)))) |
|
109 degA = degB |
|
110 |
|
111 for k in range(degB+1): |
|
112 A[i,(k*cA)+j] = B[0,k] |
|
113 |
|
114 if degA > degB: |
|
115 for k in range(degB+1,degA+1): |
|
116 A[i,(k*cA)+j] = 0 |
|
117 A, degA = clcoef(A,degA) |
|
118 return A, degA |
|
119 |
|
120 |
|
121 def ext(A, degA, k, l): |
|
122 from clcoef import clcoef |
|
123 rA, cA = polsize(A, degA) |
|
124 degB = degA |
|
125 B = pl.zeros((1, degB+1)) |
|
126 for m in range(degB+1): |
|
127 B[0, m] = A[k, (m*cA)+l] |
|
128 B,degB = clcoef(B, degB) |
|
129 return B, degB |
|
130 |
|
131 |
|
132 def transp(Q, degQ): |
|
133 """ Function to transpose a polynomial matrix. """ |
|
134 rQ, cQ = polsize(Q, degQ) |
|
135 rP = cQ |
|
136 cP = rQ |
|
137 degP = degQ |
|
138 P = pl.zeros((rP, (degP+1)*cP)) |
|
139 for i in range(degP+1): |
|
140 P[:, i*cP:(i+1)*cP] = Q[:, i*cQ:(i+1)*cQ].T |
|
141 |
|
142 return P, degP |
|
143 |
|
144 |
|
145 if __name__== "__main__": |
|
146 |
|
147 # print "Test for polsize" |
|
148 # print polsize([1, 2, 1],4) |
|
149 |
|
150 # print "Test for polmul" |
|
151 # C = pl.array([[1, 0, 0.5, 2], [0, 1, -4.71, 2.8]]) |
|
152 # A = pl.array([0.5, 3.5]) |
|
153 # print polmul(A, 0, C, 1) |
|
154 |
|
155 # print "Test for polsplit3" |
|
156 # print polsplit3([1, -0.37]) |
|
157 |
|
158 print "Test for putin" |
|
159 A = pl.array([0,0]) |
|
160 B = pl.array([0.44, -1.6, 1.6, -0.44]) |
|
161 print putin(A, 0, B, 3, 0, 0) |
|
162 |
|
163 |
|
164 pass |