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


if platform.system().upper().startswith("WIN"): 
    shared_lib_ext = ".dll"
else:
    shared_lib_ext = ".so"

# handles import from site-packages as well as local import:
_fullpath=os.path.join(os.path.dirname(__file__), "lsei"+shared_lib_ext)

# CDLL() loads the library, which can be accessed via _lib from now on
_lib = c.CDLL(_fullpath)


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

    if E != None:
        assert f != None
        assert n == E.shape[1]
    else:
        assert f == None
        E = zeros((0, n))
        f = zeros((0, 1))
    if G != None:
        assert h != None
        assert n == G.shape[1]
    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))
    pv = lambda f: f.flatten()[:, newaxis]
    Wr = vstack((pv(f), pv(b), pv(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)
    
    #activate automatic column scaling of A

    PRGOPT = array((4,    # link to next node = PRGOPT(4) using fortran indexing
                    2,    # key = 2: set/unset automatic columns scaling 
                    1,    # value = 1: set scaling
                    1),   # link = 1 --> end of list
                    dtype = float_)   

    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
              )

    #     MODE          INTEGER FLAG THAT INDICATES THE SUBPROGRAM
    #                   STATUS AFTER COMPLETION.  IF MODE.GE.2, NO
    #                   SOLUTION HAS BEEN COMPUTED.
    #
    #                   MODE =
    #
    #                   0  BOTH EQUALITY AND INEQUALITY CONSTRAINTS
    #                      ARE COMPATIBLE AND HAVE BEEN SATISFIED.
    #
    #                   1  EQUALITY CONSTRAINTS ARE CONTRADICTORY.
    #                      A GENERALIZED INVERSE SOLUTION OF EX=F WAS USED
    #                      TO MINIMIZE THE RESIDUAL VECTOR LENGTH F-EX.
    #                      IN THIS SENSE, THE SOLUTION IS STILL MEANINGFUL.
    #
    #                   2  INEQUALITY CONSTRAINTS ARE CONTRADICTORY.
    #
    #                   3  BOTH EQUALITY AND INEQUALITY CONSTRAINTS
    #                      ARE CONTRADICTORY.
    #
    #                   THE FOLLOWING INTERPRETATION OF
    #                   MODE=1,2 OR 3 MUST BE MADE.  THE
    #                   SETS CONSISTING OF ALL SOLUTIONS
    #                   OF THE EQUALITY CONSTRAINTS EX=F
    #                   AND ALL VECTORS SATISFYING GX.GE.H
    #                   HAVE NO POINTS ON COMMON.  (IN
    #                   PARTICULAR THIS DOES NOT SAY THAT
    #                   EACH INDIVIDUAL SET HAS NO POINTS
    #                   AT ALL, ALTHOUGH THIS COULD BE THE
    #                   CASE.)
    #
    return MODE.value, X


if __name__ == "__main__":

    A = eye(2) # ((1,2))
    b = zeros((2,))

    G = None
    h = None

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

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