#encoding:latin-1

# for COPYRIGHT and LICENSE see LICENSE file !

from numpy import require, zeros, double
from numpy.ctypeslib import ndpointer
from ctypes import POINTER, byref, c_int, c_double, CDLL
import os, platform

###########################################################################
#
# load library
#
# _fullpath as constructed below will find point at the libraray as long 
# as this script and the library are in the same folder. no matter if this
# folder is a local folder for development or if it is pythons public
# folder for modules aka .../Lib/site-packages/
#
###########################################################################

if platform.system().upper().startswith("WIN"): 
    _fullpath=os.path.join(os.path.dirname(__file__),"./nnls.dll")
else:
    _fullpath=os.path.join(os.path.dirname(__file__),"./nnls.so")

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

###########################################################################
#
# now we configure NNLS procedure argument types. some remarks:
#
#  *  NNLS in FORTRAN is accessed by _lib.nnls_ (the latter underscore is
#  FORTRAN convention: a function/routine which is name "xxx" in FORTRAN, 
#  will be accessed as "xxx_" in the library.
# 
#  *  "ndpointer" declares numpy array parameter, the parameters below are
#  easy to understand,  flags=_fortran_memory_layout declares FORTRAN
#  memory layout
#
#  *  parmeters declared in FORTRAN as INTEGER, DOUBLE PRECISION or FLOAT
#  are called by reference. As an example we declare an FORTRAN INTEGER
#  argument as POINTER(c_int).
#
###########################################################################

_fortran_memory_layout = [ "f_contiguous", "aligned", "writeable" ]

# look at NLLS.F for comparison:

_lib.nnls_.argtypes = [ 
    ndpointer(double, ndim=2, flags=_fortran_memory_layout),   # A      
    POINTER(c_int),                                            # MDA    
    POINTER(c_int),                                            # M      
    POINTER(c_int),                                            # N     
    ndpointer(double, ndim=1, flags=_fortran_memory_layout),   # B     
    ndpointer(double, ndim=1, flags=_fortran_memory_layout),   # X     
    POINTER(c_double),                                         # RNORM
    ndpointer(double, ndim=1, flags=_fortran_memory_layout),   # W   
    ndpointer(double, ndim=1, flags=_fortran_memory_layout),   # ZZ 
    ndpointer(int   , ndim=1, flags=_fortran_memory_layout),   # INDEX 
    POINTER(c_int),                                            # MODE 
]

###########################################################################
#
# here comes the python wrapper
#
###########################################################################

def nnls(A, b):

    """
         solves 
                  || Ax-b || -> min 

         under the constraint 

                   x >= 0 

         returns ok-flag, solution and residual:

             isok, x, res = nnls(A,b)
    """

    # 1) check, copy and convert inputs A and b to fit the declaration in 
    #    _lib.nnls_.argtypes (see above)
    M,N = A.shape
    if  M!=b.shape[0]:
        raise RuntimeError("A and b do not fit")

    A = require(A.copy(), double, _fortran_memory_layout)  # copy() because
                                                           # NNLP changes A
    B = require(b.copy(), double, _fortran_memory_layout)  # and B

    # 2) allocate working space W, ZZ, INDEX to fit declaration above #####
    W = require(zeros((N,), dtype=double), None, _fortran_memory_layout)
    ZZ= require(zeros((M,), dtype=double), None, _fortran_memory_layout)
    INDEX= require(zeros((N,), dtype=int), None, _fortran_memory_layout)

    # 3) allocate space for solution vector X  ############################
    X = require(zeros((N,), dtype=double), None, _fortran_memory_layout) 

    # 4) declare M, N, RNORM, MODE as declared above ######################
    M=c_int(M)
    N=c_int(N)
    RNORM = c_double(0.0)
    MODE=c_int()

    # 5) call NNLS  #######################################################
    _lib.nnls_(A, byref(M), byref(M), byref(N), B, X, byref(RNORM),
               W, ZZ, INDEX, byref(MODE))

    # 6) extract solution and further information and return it ###########
    return MODE.value==1, X, RNORM.value



###########################################################################
#
# and now a test ......................
#
###########################################################################
if __name__ == "__main__":

    from numpy import *
    A = random.rand(5,3)
    x = arange(3)
    b = dot(A,x)

    isok, X, res = nnls(A,b)
    
    print "X=", X
    print "res=", res
    print "residual vector= ", dot(A, X)-b