pyasl.pyx 12.8 KB
# This Source Code Form is subject to the terms of the Mozilla Public
# License, v. 2.0. If a copy of the MPL was not distributed with this
# file, You can obtain one at http://mozilla.org/MPL/2.0/.

cimport casl, cpython, libc.stdio
from libc.stdlib cimport malloc, calloc, free
import numpy as np
import os, scipy.sparse

cdef enum:
    ASL_read_f    = 1  # no derivatives
    ASL_read_fg   = 2  # first derivatives
    ASL_read_fgh  = 3  # first derivatives and hessian vector products
    ASL_read_pfg  = 4  # first derivatives and partially seperable structure
    ASL_read_pfgh = 5  # first and second derivatives and partially seperable structure


cdef class ASL:
    cdef casl.ASL* asl
    cdef libc.stdio.FILE* nl
    # Table 1 from Hooking your solver...
    cdef int asl_n_var              # number of variables
    cdef int asl_nbv                # number of binary variables
    cdef int asl_niv                # number of integer variables
    cdef int asl_n_con              # number of constraints
    cdef int asl_n_obj              # number of objectives
    cdef int asl_nlo                # number of nonlinear objectives, appear before linear
    cdef int asl_nranges            # number of ranged constraints (finite boxes)
    cdef int asl_nlc                # number of nonlinear general constraints
    cdef int asl_nlnc               # number of nonlinear network constraints
    cdef int asl_nlvb               # number of variables appearing nonlinearly in constraints and objectives
    cdef int asl_nlvbi              # number of integer variables appearing nonlinearly in constraints and objectives
    cdef int asl_nlvc               # number of variables appearing nonlinearly in constraints
    cdef int asl_nlvci              # number of integer variables appearing nonlinearly in constaints
    cdef int asl_nlvo               # number of variables appearing nonlinearly in objectives
    cdef int asl_nlvoi              # number of integer variables appearing nonlinearly in objectives
    cdef int asl_lnc                # number of linear network constraints
    cdef int asl_nzc                # number of nonzeros in jacobian
    cdef int asl_nzo                # number of nonzeros in objective gradients
    cdef int asl_maxrownamelen      # length of longest constraint or objective name
    cdef int asl_maxcolnamelen      # length of longest variable name
    # Table 2
    cdef double* asl_x0             # primal initial guess
    cdef double* asl_LUv            # array of lower (and, if UVx is null, upper) variable bounds
    cdef double* asl_Uvx            # array of upper variable bounds
    cdef double* asl_LUrhs          # array of lower (and, if Urhsx is null, upper) constraint bounds
    cdef double* asl_Urhsx          # array of upper constraint bounds
    # complementarity information
    cdef int asl_n_cc                # number of complementarity constraints
    cdef int* asl_cvar              # cvar[i] > 0 means constraint i complements variable cvar[i] - 1
    
    cdef int* colind      
    cdef int* rowind      
    
    def __cinit__(self):
        self.asl = casl.ASL_alloc(ASL_read_pfgh)
        if self.asl is NULL:
            cpython.PyErr_NoMemory()
            
    def __init__(self, stub):
        # check that stub file really exists
        base, ext = os.path.splitext(stub)
        if ext.lower() != '.nl':
            stub += '.nl'
        
        if not os.path.exists(stub):
            raise Exception("File %s does not exist!" % stub)
        
        self.asl.i.want_deriv_ = 1;
        self.asl.i.return_nofile_ = 1;
        self.asl.i.want_xpi0_ = 1;        
        self.asl.i.Fortran_ = 0;

        # give stub to ASL
        self.nl = casl.jac0dim_ASL(self.asl, stub, len(stub))
        
        # read quantities from Table 1
        self.asl_n_var         = self.asl.i.n_var_
        self.asl_nbv           = self.asl.i.nbv_          
        self.asl_niv           = self.asl.i.niv_          
        self.asl_n_con         = self.asl.i.n_con_        
        self.asl_n_obj         = self.asl.i.n_obj_        
        self.asl_nlo           = self.asl.i.nlo_          
        self.asl_nranges       = self.asl.i.nranges_ 
        self.asl_nlc           = self.asl.i.nlc_         
        self.asl_nlnc          = self.asl.i.nlnc_        
        self.asl_nlvb          = self.asl.i.nlvb_        
        self.asl_nlvbi         = self.asl.i.nlvbi_       
        self.asl_nlvc          = self.asl.i.nlvc_        
        self.asl_nlvci         = self.asl.i.nlvci_       
        self.asl_nlvo          = self.asl.i.nlvo_        
        self.asl_nlvoi         = self.asl.i.nlvoi_       
        self.asl_lnc           = self.asl.i.lnc_         
        self.asl_nzc           = self.asl.i.nzc_         
        self.asl_nzo           = self.asl.i.nzo_         
        self.asl_n_cc          = self.asl.i.n_cc_
        self.asl_maxrownamelen = self.asl.i.maxrownamelen_
        self.asl_maxcolnamelen = self.asl.i.maxcolnamelen_
        
        # if self.asl_n_obj > 1:
        #     raise Exception("Multiobjective Optimization not supported")
                
        # allocate storage for initial starting point and bounds
        self.asl_x0       = <double *> malloc( self.asl_n_var * sizeof(double) )
        self.asl.i.X0_    = self.asl_x0
        self.asl_LUv      = <double *> calloc( self.asl_n_var, sizeof(double) )
        self.asl.i.LUv_   = self.asl_LUv
        self.asl_Uvx      = <double *> calloc( self.asl_n_var, sizeof(double) )
        self.asl.i.Uvx_   = self.asl_Uvx
        self.asl_LUrhs    = <double *> malloc( self.asl_n_con * sizeof(double) )
        self.asl.i.LUrhs_ = self.asl_LUrhs
        self.asl_Urhsx    = <double *> malloc( self.asl_n_con * sizeof(double) )
        self.asl.i.Urhsx_ = self.asl_Urhsx
        
        # allocate storage for complementarity mapping
        if self.asl_n_cc > 0:
            self.asl_cvar = <int *> malloc( self.asl_n_con * sizeof(int) )
            self.asl.i.cvar_ = self.asl_cvar

        # get the full nonlinear problem
        flags = 12 | 16 | 0x80 # ASL_findgroups | ASL_return_read_err | ASL_keep_all_suffixes
        casl.pfgh_read_ASL(self.asl, self.nl, flags)
        
        # if there are constraints get sparsity structure of jacobian
        cdef int ii
        if self.asl_n_con > 0:
            self.colind  = <int *> malloc( self.asl_nzc * sizeof(int) )
            self.rowind  = <int *> malloc( self.asl_nzc * sizeof(int) )
            for ii in range(self.asl_n_con):
                cgrad = self.asl.i.Cgrad_[ii]
                while cgrad is not NULL:
                    self.rowind[cgrad.goff] = ii
                    self.colind[cgrad.goff] = cgrad.varno
                    cgrad = cgrad.next
        else:
            self.colind = NULL
            self.rowind = NULL
        
    property n:
        """number of variables"""
        def __get__(self):
            return self.asl_n_var
        
    property m:
        """number of constraints"""
        def __get__(self):
            return self.asl_n_con
        
    property ncc:
        def __get__(self):
            return self.asl_n_cc

    property x0:
        """starting point"""
        def __get__(self):
            cdef double [:] tmp = <double[:self.asl_n_var]> self.asl_x0
            return np.asarray(tmp)
        
    property xl:
        """lower variable bounds"""
        def __get__(self):
            cdef double [:] tmp = <double[:self.asl_n_var]> self.asl_LUv
            return np.asarray(tmp)
        
    property xu:
        """upper variable bounds"""
        def __get__(self):
            cdef double [:] tmp = <double[:self.asl_n_var]> self.asl_Uvx
            return np.asarray(tmp)

    property cl:
        """lower constraints bounds"""
        def __get__(self):
            cdef double [:] tmp
            if self.asl_n_con > 0:
                tmp = <double[:self.asl_n_con]> self.asl_LUrhs
                return np.asarray(tmp)
            return np.array([])
        
    property cu:
        """upper constraints bounds"""
        def __get__(self):
            cdef double [:] tmp
            if self.asl_n_con > 0:
                tmp = <double[:self.asl_n_con]> self.asl_Urhsx
                return np.asarray(tmp)
            return np.array([])
        
    property cvar:
        """complementarity mapping"""
        def __get__(self):
            cdef int [:] tmp
            if self.asl_n_cc > 0:
                tmp = <int[:self.asl_n_con]> self.asl_cvar
                return np.asarray(tmp)
            return np.array([])
        
    def obj(self, double[:] x_view):
        """Objective value at x"""
        if self.asl_n_obj == 0:
            return 0.
        cdef casl.fint error = 0
        cdef double objval
        objval = (self.asl.p.Objval)(self.asl, 0, &x_view[0], &error)
        if error != 0:
            raise Exception("AMPL could not evaluate objective")
        else:
            return objval
        
    def con(self, double[:] x_view):
        """Constraints vector at x"""
        cdef casl.fint error = 0
        c = np.zeros([self.asl_n_con])
        if self.asl_n_con == 0:
            return c
        cdef double [:] c_view = c
        (self.asl.p.Conval)(self.asl, &x_view[0], &c_view[0], &error)
        if error != 0:
            raise Exception("AMPL could not evaluate constraints")
        else:
            return c
        
    def grad(self, double[:] x_view):
        """Gradient of Objective at x"""
        if self.asl_n_obj == 0:
            return 0.
        cdef casl.fint error = 0
        grad = np.zeros([self.asl_n_var])
        if self.asl_n_obj == 0:
            grad[:] = 0
            return grad
        cdef double [:] grad_view = grad
        (self.asl.p.Objgrd)(self.asl, 0, &x_view[0], &grad_view[0], &error)
        if error != 0:
            raise Exception("AMPL could not evaluate gradient")
        else:
            return grad
        
    def jac(self, double[:] x_view):
        """Jacobian of Constraints at x"""
        if self.asl_n_con == 0:
            return scipy.sparse.coo_matrix((self.asl_n_con, self.asl_n_var))
        jacdata = np.zeros([self.asl_nzc])
        cdef double [:] jacdata_view = jacdata
        cdef casl.fint error
        (self.asl.p.Jacval)(self.asl, &x_view[0], &jacdata_view[0], &error)
        cdef int[:] colview = <int[:self.asl_nzc]> self.colind
        cdef int[:] rowview = <int[:self.asl_nzc]> self.rowind
        if error != 0:
            raise Exception("AMPL could not evaluate jacobian")
        else: 
            return scipy.sparse.coo_matrix((jacdata, (np.asarray(rowview), np.asarray(colview) )), shape=(self.asl_n_con, self.asl_n_var))
        
    def hessprod(self, l, sigma, d):
        """Computes [sigma*f'' + sum_i l_i * c_i'']*d at the x where objective, constraints or their derivative have been computed"""
        cdef double* ow = NULL
        if self.asl_n_obj == 1:
            ow = <double*> malloc( sizeof(double) )
            ow[0] = sigma
        hessprod = np.zeros([self.asl_n_var])
        cdef double[:] hessprod_view = hessprod
        cdef double[:] l_view = l
        cdef double[:] d_view = d
        (self.asl.p.Hvcomp)(self.asl, &hessprod_view[0], &d_view[0], 0, &ow[0], &l_view[0] if self.asl_n_con > 0 else NULL)
        if self.asl_n_obj == 1:
            free(ow)
        return hessprod

    def fulhes(self, int no, double ow, double[:] y_view):
        """Computes full hessian [ow*f'' + sum_i y_i * c_i''] at x where objective, constraints or their derivative have been computed"""
        
        hess = np.zeros([self.asl_n_var,self.asl_n_var])
        cdef double[:,:] hess_view = hess
        cdef double ow_view = ow
        (self.asl.p.Fulhes)(self.asl, &hess_view[0,0], self.asl_n_var, no, &ow_view, &y_view[0] if self.asl_n_con > 0 else NULL)
        return hess

    def writesol(self, msg = '', x=None, y=None):
        """writes ampl .sol file, gives back termination message msg"""
        cdef double [:] x_view
        cdef double [:] y_view
        if x != None:
            assert x.shape[0] == self.asl_n_var
            x_view = x
        if y != None:
            assert y.shape[0] == self.asl_n_con
            y_view = y
        casl.write_sol_ASL(self.asl, msg, &x_view[0] if x != None else NULL, &y_view[0] if y != None else NULL, NULL)
        
         
    def __dealloc__(self):
        return
        # free(self.asl_x0)
        # free(self.asl_LUv)
        # free(self.asl_Uvx)
        # free(self.asl_LUrhs)
        # free(self.asl_Urhsx)
        # if self.colind is not NULL:
        #     free(self.colind)
        # if self.rowind is not NULL:
        #     free(self.rowind)
        # if self.asl is not NULL:
        #     free(self.asl)