|
1 #!/usr/bin/python |
|
2 |
|
3 import pylab as pl |
|
4 from rowjoin import rowjoin |
|
5 from polyfuncs import polsize |
|
6 from t1calc import t1calc |
|
7 from seshft import seshft |
|
8 from colsplit import colsplit |
|
9 from clcoef import clcoef |
|
10 |
|
11 def left_prm(N, degN, D, degD, job=1, gap=None): |
|
12 """ """ |
|
13 if gap == None: |
|
14 gap = 1.0e8 |
|
15 |
|
16 F, degF = rowjoin(D, degD, N, degN) |
|
17 Frows, Fbcols = polsize(F, degF) |
|
18 Fcols = Fbcols * (degF+1) |
|
19 T1 = pl.array([]) |
|
20 pr = pl.array([]) |
|
21 degT1, T1rows, shft = 0, 0, 0 |
|
22 S=F |
|
23 sel = pl.ones((Frows,1)) |
|
24 T1bcols =1 |
|
25 abar = pl.arange(Fbcols, Frows) |
|
26 if T1rows: |
|
27 T1 = pl.hstack((T1, pl.zeros((T1rows, Frows)))) |
|
28 if T1.size == 0: |
|
29 T1rows = 0 |
|
30 while T1.size == 0 or T1rows < (Frows - Fbcols): |
|
31 if T1.size == 0: |
|
32 T1rows = 0 |
|
33 Srows = Frows*T1bcols |
|
34 T1, T1rows, sel, pr = t1calc(S, Srows, T1, T1rows, sel, pr, Frows, Fbcols, abar, gap) |
|
35 T1rows, T1cols = pl.atleast_2d(T1).shape |
|
36 if T1.size == 0: |
|
37 T1rows = 0 |
|
38 if T1rows < (Frows - Fbcols): |
|
39 if T1.size != 0: |
|
40 T1 = pl.hstack((T1, pl.zeros((T1rows,Frows)).squeeze())) |
|
41 T1bcols += 1 |
|
42 degT1 += 1 |
|
43 shft += Fbcols |
|
44 S = seshft(S,F,shft); |
|
45 sel = pl.hstack((sel,sel[Srows-Frows:Srows+1])) |
|
46 rowvec = pl.arange((T1bcols-1)*Frows+(Fbcols), (T1bcols * Frows)) |
|
47 abar = pl.hstack((abar, rowvec)) |
|
48 |
|
49 B, degB, A, degA = colsplit(T1, degT1, Fbcols, Frows-Fbcols) |
|
50 B, degB = clcoef(B, degB) |
|
51 B = -B; |
|
52 A, degA = clcoef(A, degA) |
|
53 if job == 1: |
|
54 return B, degB, A, degA |
|
55 |
|
56 if job == 2: |
|
57 S = S[sel!=0,:] |
|
58 redSrows, Scols = S.shape |
|
59 C = pl.hstack(pl.eye(Fbcols), pl.zeros(Fbcols,Scols-Fbcols)) |
|
60 T2 = pl.dot(C, pl.inv(S)) |
|
61 T2 = makezero(T2,gap); |
|
62 T2 = move(T2, pl.find(sel), Srows); |
|
63 X, degX, Y, degY = colsplit(T2, degT1, Fbcols, Frows - Fbcols); |
|
64 X, degX = clcoef(X, degX); |
|
65 Y, degY = clcoef(Y, degY); |
|
66 elif job == 3: |
|
67 Y = S |
|
68 degY = sel |
|
69 X = degT1 |
|
70 degX = Fbcols |
|
71 |
|
72 else: |
|
73 print 'Message from left_prm:no legal job number specified' |
|
74 exit() |
|
75 |
|
76 return B, degB, A, degA, Y, degY, X, degX |
|
77 |
|
78 |
|
79 if __name__== "__main__": |
|
80 D = pl.array([[1, 0, 0, 0, 0, 0], |
|
81 [0, 1, 0, 1, 0, 0], |
|
82 [0, 0, 1, 1, 1, 0]]) |
|
83 |
|
84 N = pl.array([[1, 0, 0], |
|
85 [0, 1, 0], |
|
86 [0, 0, 1]]) |
|
87 |
|
88 dD, dN = 1, 0 |
|
89 |
|
90 result = left_prm(N, dN, D, dD) |
|
91 print result |
|
92 pass |