#encoding:latin-1
from   numpy import *
import ctypes as c

lib = c.cdll["lsei.dll"]

def LSEI(A, b, E=None, f=None, G=None, h=None):
    n = A.shape[1]

    pv = lambda f: f.flatten()[:, newaxis]
    b = pv(b)

    assert b.shape[0] == A.shape[0]

    if E != None:
        assert f != None
        assert n == E.shape[1]
        f = pv(f)
        assert f.shape[0] == E.shape[0]
    else:
        assert f == None
        E = zeros((0, n))
        f = zeros((0, 1))

    if G != None:
        assert h != None
        assert n == G.shape[1]
        h = pv(h)
        assert h.shape[0] == G.shape[0]
    else:
        assert h == None
        G = zeros((0, n))
        h = zeros((0, 1))

    ma = A.shape[0]
    me = E.shape[0]
    mg = G.shape[0]
    WM = vstack((E, A, G))
    Wr = vstack((f, b, h))
    W =  hstack((WM, Wr))

    lip = mg+2*n+2
    ip = zeros((lip,), dtype=int_)
    k = max(ma+mg, n)
    lw = 2*(me+n)+k+(mg+2)*(n+7)
    ip[0] = lw
    ip[1] = lip

    req = "F_CONTIGUOUS ALIGNED".split()
    W = require(W, float_, req)

    MDW = c.c_int(ma+me+mg)
    ME  = c.c_int(me)
    MA  = c.c_int(ma)
    MG  = c.c_int(mg)
    N   = c.c_int(n)
    PRGOPT = array((1.0,))
    X = zeros((n,))
    RNORME = c.c_double()
    RNORML = c.c_double()
    MODE   = c.c_int()
    WS     = zeros((lw, ))
    IP     = zeros((lip, ), dtype=int_)

    lib.lsei_.argtypes = [ \
                            ctypeslib.ndpointer(float_, ndim=2, flags=req),  # W
                            c.POINTER(c.c_int),   # MDW
                            c.POINTER(c.c_int),   # ME
                            c.POINTER(c.c_int),   # MA
                            c.POINTER(c.c_int),   # MG
                            c.POINTER(c.c_int),   # N
                            ctypeslib.ndpointer(float_, ndim=1), # PROGOPT
                            ctypeslib.ndpointer(float_, ndim=1, flags="contiguous, writeable"), # X
                            c.POINTER(c.c_double), # RNORME
                            c.POINTER(c.c_double), # RNORML
                            c.POINTER(c.c_int),    # MODE
                            ctypeslib.ndpointer(float_, ndim=1, flags="contiguous"), # WS
                            ctypeslib.ndpointer(int_, ndim=1, flags="contiguous"), # IP
                         ]

    lib.lsei_(  W, 
                c.byref(MDW), 
                c.byref(ME),
                c.byref(MA),
                c.byref(MG),
                c.byref(N),
                PRGOPT,
                X,
                c.byref(RNORME),
                c.byref(RNORML),
                c.byref(MODE),
                WS, 
                IP
             )

    return MODE.value, X

if __name__ == "__main__":

    n = 12
    m = 20

    A = random.rand(m, n)
    b = zeros((m,))

    print linalg.svd(A)[1]

    G = -eye(n)
    h = -0.5 * ones((n,))

    E = ones((1,n))
    f = ones((1,))

    mode, x = LSEI(A, b, E, f, G, h)
    print "mode = ", mode
    print "x = ", x
    print sum(x)