python/xdync.py
changeset 0 0efde00f9229
equal deleted inserted replaced
-1:000000000000 0:0efde00f9229
       
     1 #!/usr/bin/python
       
     2 
       
     3 import pylab as pl
       
     4 from rowjoin import rowjoin
       
     5 from polyfuncs import polsize
       
     6 from cindep import cindep
       
     7 from seshft import seshft
       
     8 from colsplit import colsplit
       
     9 from left_prm import left_prm
       
    10 from move import move
       
    11 from clcoef import clcoef
       
    12 
       
    13 def xdync(N, degN, D, degD, C, degC, gap=1.0e8):
       
    14     """ """
       
    15     C = pl.atleast_2d(C)
       
    16     F, degF = rowjoin(D, degD, N, degN)
       
    17     Frows, Fbcols = polsize(F, degF)
       
    18     B, degB, A, degA, S, sel, degT1, Fbcols = left_prm(N, degN, D, degD, 3, gap)
       
    19     Crows, Ccols = C.shape
       
    20     Srows, Scols = pl.atleast_2d(S).shape
       
    21     S = S[sel!=0,:]
       
    22     T2 = pl.array([])
       
    23 
       
    24     for i in range(Crows):
       
    25         Saug = seshft(S, C[i,:], 0)
       
    26         b = cindep(Saug)
       
    27         b = move(b, pl.find(sel), Srows)
       
    28         if T2.size != 0:
       
    29             T2 = pl.vstack((T2, b))
       
    30         else:
       
    31             T2 = b.copy()
       
    32 
       
    33     X, degX, Y, degY = colsplit(T2, degT1, Fbcols, Frows-Fbcols)
       
    34     X, degX = clcoef(X, degX)
       
    35     Y, degY = clcoef(Y, degY)
       
    36 
       
    37     return Y, degY, X, degX, B, degB, A, degA
       
    38 
       
    39 
       
    40 if __name__== "__main__":
       
    41     N = pl.array([[0,  4, 0, 1],
       
    42                [-1, 8, 0, 3]])
       
    43     dN = 1
       
    44     D = pl.array([[0, 0,  1, 4, 0, 1],
       
    45                [0, 0, -1, 0, 0, 0]])
       
    46     dD = 2
       
    47     C = pl.array([[1, 0, 1, 1],
       
    48                [0, 2, 0, 1]])
       
    49     dC = 1
       
    50 
       
    51     print xdync(N,dN,D,dD,C,dC)
       
    52 
       
    53