0
|
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
|